Módulo 3 - Dados de Área

Neste módulo vamos mostrar as principais técnicas de análise estatística espacial utilizando dados de área na vigilância em saúde.

Ao final deste módulo, você será capaz de:

  1. Importar e manipular dados de área no software R.
  2. Construir mapas atentando-se a suas escalas, pontes de corte e legendas.
  3. Entender o conceito de autocorrelação espacial, assim como identificá-la global e localmente dentro de um conjunto de dados espacial.

Criando mapas

A construção de mapas com dados de área é muito útil na vigilância em saúde pública pois permite visualizar espacialmente a distribuição em áreas de doenças, recursos, fatores de risco, entre outros. Alguns cuidados são necessários para que tais mapas não sejam enviesados, ou seja, que a sua visualização não esteja induzindo a interpretações incorretas. A seguir, vamos mostrar como a definição dos pontos de corte de uma variável contínua (por exemplo, incidência de uma doença) pode levar a conclusões enviesadas, e quais são as boas práticas para evitá-las.

Pontos de corte em mapas

Ao representar uma variável em mapa (exemplo: número de casos, ou incidência de casos de uma doença), nem sempre representá-la de forma contínua pode ser informativo. Em distribuções desbalanceadas e/ou que possuem valores aberrantes (valores próximos de zero e valores muito mais altos que os demais), a escala de cores para o mapa pode acabar chamando atenção somente para valores extremos.

Vamos tomar como exemplo a incidência de casos prováveis de dengue por município no estado de São Paulo no ano de 2015. A partir do número de casos notificados no SINAN e da população estimada para o ano em questão, podemos calcular a incidência e representá-la, de forma contínua, em um mapa como na Figura 18.

Figura 18: Incidência de casos prováveis de dengue por município no estado de São Paulo, 2015.

Figura 18: Incidência de casos prováveis de dengue por município no estado de São Paulo, 2015.

No mapa, nota-se alguns municípios que chamam mais atenção pela cor mais escura - ou seja, maior incidência. Alguns outros pontos de concentração da doença são marcados no estado, mas, no geral, a maioria dos municípios é pintada em tom mais claro.

A visualização da variável em sua forma contínua acaba sendo útil em poucos casos, principalmente em situações ligadas à vigilância em saúde. Tradicionalmente, trabalhamos com dados de contagem de casos - e, na maioria dos cenários, a distribuição dessa variável é desbalanceada. Poucos municípios possuem muitos casos enquanto muitos municípios possuem poucos casos. Além disso, municípios com populações pequenas acabam tendo a incidência inflada por conta de um denominador pequeno. Isso faz com que poucos municípios se destaquem dos demais na escala de cores utilizando a variável de forma contínua, e torna difícil a distinção do cenário entre o restante.

Dessa forma, torna-se importante a separação dessa variável em faixas - através da definição de pontos de corte. Contudo, a escolha da estratégia utilizada para definição desses pontos importa e pode impactar e muito no resultado da visualização obtida.

Definindo pontos de corte

Observe os mapas abaixo na Figura 19. Eles representam a mesma informação - incidência de casos prováveis de dengue por município no estado de São Paulo, em 2015. O que há, portanto, de diferente entre eles?

Figura 19: Formas de representação da incidência de casos prováveis de dengue no estado de São Paulo, 2015.

Figura 19: Formas de representação da incidência de casos prováveis de dengue no estado de São Paulo, 2015.

Podemos observar que, a depender da estratégia de definição dos pontos de corte utilizada, a interpretação que temos de cada mapa pode ser diferente. Em alguns mapas (A e C), a situação no estado parece pior, com mais municípios pintados em tons escuros, levando a pensar em uma maior incidência. No mapa B, no entanto, a impressão é exatamente a oposta - a maioria dos municípios está em tons claros, o que nos leva a interpretar que o estado em geral está com baixa incidência da doença. Vamos entender cada uma dessas situações.

Quebras arbitrárias

O primeiro caso, representado na Figura 19A, não utiliza nenhum método específico para estabelecimento de pontos de corte da incidência. Os pontos são definidos arbitrariamente, e podem ou não ser baseados em uma referência existente sobre o tema. Em um departamento de vigilância, pode-se ter estabelecido que uma incidência de 50 casos por 10.000 habitantes é um patamar controlado, e acima de 100 casos por 100.000 tem-se um cenário de preocupação. Assim, o analista pode-se basear nesses pontos pré-definidos na rotina de vigilância para gerar seus mapas. Contudo, não é garantia de que a visualização seja adequada aos dados: pode ser que muitos municípios se concentrem em uma mesma faixa, por exemplo, trazendo dificuldades na distinção entre esses municípios (Figura 20).

Figura 20: Distribuição dos municípios em cada faixa de incidência, definida por quebras arbitrárias.

Figura 20: Distribuição dos municípios em cada faixa de incidência, definida por quebras arbitrárias.

Quebras regulares

A segunda estratégia de quebras utilizada é a das quebras regulares. Nela, considera-se o intervalo entre o mínimo e o máximo da distribuição, e divide-se esse intervalo em partes iguais. Dessa forma, cada faixa das quebras possui a mesma amplitude.

Repare, na Figura 19B, que a distância entre um ponto de corte e outro é de aproximadamente 296 para todos os intervalos.

Essa estratégia facilita a comunicação, mas é particularmente útil quando a distribuição dos dados se aproxima de uma distribuição uniforme. No cenário da vigilância em saúde, esse é raramente o caso - há muitos valores concentrados em uma faixa da distribuição e poucos valores distribuídos entre outras. Isso causa o efeito de que um número excessivo de municípios pode cair na mesma categoria, enquanto outras categorias podem possuir um ou dois municípios (Figura 21). Pode ser, até, que faixas de incidência não possuam nenhum município dentro dela.

Figura 21: Distribuição dos municípios em cada faixa de incidência, definida por quebras regulares.

Figura 21: Distribuição dos municípios em cada faixa de incidência, definida por quebras regulares.

Vê-se agora que a distância entre a mudança de cores é sempre a mesma, o que é agradável visualmente. Contudo, a vasta maioria dos municípios se concentra dentro de uma mesma cor. Dessa forma, a estratégia falha em distinguir esses municípios com incidência mais baixa entre eles.

Quebras quantílicas

Outra estratégia consiste em determinar pontos de quebra a partir de quantis da distribuição. As quebras quantílicas tentam dividir os municípios em grupos de forma com que cada grupo contenha aproximadamente a mesma quantidade de municípios. Por exemplo, ao dividir em 5 grupos, teríamos cada grupo com aproximadamente 20% do total de observações.

Figura 22: Distribuição dos municípios em cada faixa de incidência, definida por quebras quantílicas

Figura 22: Distribuição dos municípios em cada faixa de incidência, definida por quebras quantílicas

Dessa forma, “forçamos” os grupos a conterem aproximadamente a mesma quantidade de municípios. Essa técnica garante que os mapas terão representação visível de todos os grupos; contudo, vemos que mesmo municípios com incidências baixas (Figura 22, primeira barra) já estão sendo categorizados para o segundo grupo, enquanto o ponto de corte para o grupo de incidência mais alta fica muito mais baixo. Isso pode levar a falsas impressões nos mapas, já que fazemos com que necessariamente uma parcela igual esteja contida nos níveis mais altos e mais baixos de incidência.

Quebras naturais, ou quebras de Jenks

O método de classificação por quebras naturais (ou quebras de Jenks) parte da ideia de obter pontos de corte que tentam, ao mãximo, minimizar a variação intra-grupo e maximizar a variância entre os grupos. Ou seja, busca-se dividir o conjunto de dados em grupos semelhantes entre si e diferentes dos demais.

Figura 23: Distribuição dos municípios em cada faixa de incidência, definida por quebras naturais de Jenks.

Figura 23: Distribuição dos municípios em cada faixa de incidência, definida por quebras naturais de Jenks.

Para este exemplo, as quebras naturais classificaram os municípios de forma equilibrada em relação às demais (Figura 23) - não houve uma hiper concentração dos municípios em uma só classe (como as quebras regulares), nem a definição de quebras que abrangiam muitos municípios diferentes dentro de si (como nas quebras quantílicas). Percebe-se que o primeiro grupo consta majoritariamente de municípios com baixa incidência (primeira barra no histograma), e o grupo mais alto inclui poucos municípios, mas que realmente apresentaram valores de incidência mais extremos na distribuição.

Comparando mapas no tempo

Outro ponto importante ao definir escalas e pontos de corte em mapas é a padronização dessas escalas ao comparar mapas em diferentes pontos no tempo.

Por exemplo, suponha que queremos comparar a incidência de dengue nos municípios do estado de São Paulo no período 2014-2017, utilizando as quebras quantílicas (Figura 24).

Figura 24: Incidência de dengue nos municípios São Paulo, ao longo de 2014-2017.

Figura 24: Incidência de dengue nos municípios São Paulo, ao longo de 2014-2017.

Os mapas gerados nos levam a interpretações sobre as áreas onde a doença se concentrou em cada ano. Podemos dizer que há certa similaridade entre os padrões espaciais, onde certas regiões do Norte e Noroeste paulista aparecendo repetidamente nas faixas de maior incidência. A impressão geral que temos é que, a distribuição da doença se deu de forma parecida ao longo dos anos. Contudo, vemos que os pontos de corte em cada mapa são criticamente diferentes: enquanto o tom mais escuro no mapa relativo à 2017 representa uma incidência maior que 5,67, o mesmo tom no mapa de 2015 representa incidência superior a 393,97.

Assim, ao comparar mapas no tempo, faz-se necessária a utilização de uma escala única para permitir a comparação das taxas não somente no espaço, mas entre os anos. Veja:

Figura 25: Incidência de dengue nos municípios São Paulo, ao longo de 2014-2017, utilizando escalas de cor independentes a cada ano.

Figura 25: Incidência de dengue nos municípios São Paulo, ao longo de 2014-2017, utilizando escalas de cor independentes a cada ano.

Agora, utilizamos a mesma técnica de classificação a partir de quebras quantílicas, mas para todo o período de estudo para determinação das quebras. Isso permite a comparação dos cenários ao longo do tempo: além de identificarmos certos padrões espaciais, vemos que o ano de 2015 foi muito mais crítico em termos de incidência de dengue do que os demais anos - principalmente 2017, quando a maioria dos municípios ficou no primeiro grupo da escala de cores.

Problemas conhecidos ao lidar com dados de área

Ao trabalhar com dados agregados a nível de área, estamos sujeitos à dependência das definições e delimitações impostas na divisão do território, sobre as quais muitas vezes não temos controle. Aqui tratamos de alguns efeitos conhecidos que podem acontecer ao lidar com esse tipo de dado, e que devem sempre ser considerados ao interpretar os resultados de análises com dados de área.

Tamanho heterogêneo das áreas

Como mencionado, as delimitações no território frequentemente são derivadas de divisões políticas e administrativas. Entre outros efeitos, isso pode levar à geração de áreas pequenas em extensão territorial, porém com alta densidade populacional; enquanto áreas extensas concentram pouca população. Isso pode afetar nossas visualizações: as áreas maiores, além de chamarem maior atenção aos nossos olhos, estão sujeitas a uma inflação das taxas devido ao baixo denominador (população). Vamos conferir, por exemplo, a análise da taxa de alfabetização dos municípios brasileiros, segundo o Censo 2022 (Figura 26).

As áreas de principal destaque no mapa são áreas consideravelmente grandes, e que apresentam piores indicadores.

Figura 26: Taxa de Alfabetização (%) nos municípios brasileiros, segundo o Censo 2022.

Figura 26: Taxa de Alfabetização (%) nos municípios brasileiros, segundo o Censo 2022.

Vamos comparar esse padrão espacial com a densidade populacional nesses municípios (Figura 27).

Figura 27: Densidade populacional (habitantes/km²) nos municípios brasileiros, segundo o Censo 2022.

Figura 27: Densidade populacional (habitantes/km²) nos municípios brasileiros, segundo o Censo 2022.

Vemos que além de grandes, as áreas de destaque no primeiro mapa possuem baixa densidade populacional. Como a população é frequentemente o denominador nos indicadores que analisamos, quando a população é pequena essas taxas tendem a ser instáveis. Portanto, municípios de população pequena (e que, como vimos, podem ser grandes em termos de extensão territorial) estão sujeitos a uma elevação artificial nas taxas devido a essa baixa população, e acabam mascarando demais padrões espaciais quando olhamos o mapa completo.

Problema da área modificável (MAUP)

Outro problema comum é que as taxas e indicadores dependem da forma com que as áreas são delimitadas e agregadas. E se, de repente, essas delimitações são alteradas, resultados diferentes podem ser obtidos. Esse problema é conhecido como Problema da Área Modificável (Modifiable Areal Unit Problem - MAUP).

Vamos tentar ilustrar esse efeito a partir de dados hipotéticos de casos de uma doença sobre um conjunto de setores censitários (Figura 28).

Figura 28: Número de casos de uma doença hipotética sobre um conjunto de setores censitários.

Figura 28: Número de casos de uma doença hipotética sobre um conjunto de setores censitários.

Agora, vamos supor que esses dados são agregados para uma unidade espacial maior (como bairro) para serem disponibilizados. A forma com que essas delimitações de bairro são definidas a partir dos setores censitários pode mudar a visualização da incidência da doença que obtemos. Vamos comparar duas delimitações distintas na Figura 29:

Figura 29: Comparação da incidência da doença hipotética em duas divisões de território diferentes (A e B) a partir dos setores censitários.

Figura 29: Comparação da incidência da doença hipotética em duas divisões de território diferentes (A e B) a partir dos setores censitários.

Vemos que as duas divisões, embora passem a mesma informação de distribuição da doença, mostram impressões ligeiramente diferentes de onde ela se concentra. Além disso, taxas de incidência consideravelmente mais altas são obtidas a depender do zoneamento. Como o processo de agregação pode se dar em vários níveis a partir dos dados referenciados, deve-se ter atenção para esses possíveis efeitos de atenuação ou de amplificação das taxas de incidência a depender dos níveis e divisões no processo de agregação.

Gerrymandering

A variação de proporções e taxas de acordo com a mudança dos limites territoriais pode ser usado intencionalmente e com a finalidade de alterar resultados importantes - como é feito no contexto eleitoral estadunidense. Em tal país, parte do processo eleitoral envolve contar qual partido é ganhador em cada distrito (voto distrital) e a partir disso, contar quantos distritos cada partido “tem”. A forma com que os distritos são delimitados é extremamente importante, pois pode alterar os resultados obtidos no final das contas.

Observe o exemplo a seguir na Figura 30, onde temos 50 pessoas e 60% delas votaram no partido laranja, e os 40% restantes no partido roxo.

Figura 30: Exemplo de Gerrymandering - A depender da divisão dos distritos, o resultado da contagem é alterado.
Figura 30: Exemplo de Gerrymandering - A depender da divisão dos distritos, o resultado da contagem é alterado.

Perceba que, dependendo da forma com que são traçados os distritos, temos uma mudança significativa na contagem de distritos que cada partido possui: no terceiro cenário, temos todos os distritos com vitória do partido laranja. Por outro lado, mesmo com menor proporção total de votos no total da população, a quarta configuração faz com que o partido roxo tenha mais distritos - 3 contra 2 do laranja. Esse procedimento de calcular as divisões distritais com a finalidade de alterar os resultados obtidos no processo eleitoral ocorre nos Estados Unidos, e é chamado de Gerrymandering.

Prática em R: Criação de mapas temáticos

Vamos realizar nossas primeiras visualizações, utilizando os dados de casos de dengue no município de São Paulo em 2015. Para manipular dados espaciais, vamos utilizar a biblioteca sf. É importante se atentar aos pré-requisitos de instalação da biblioteca e suas dependências. Mais detalhes sobre sua instalação podem ser encontrados aqui.

library(tidyverse)
library(sf)
  • library(tidyverse): carrega o conjunto de pacotes do R chamado tidyverse, que inclui ferramentas para manipular dados, fazer gráficos e análises estatísticas de forma prática (ex.: dplyr, ggplot2, readr etc.).

  • library(sf): carrega o pacote sf, que é usado para trabalhar com dados espaciais (mapas, coordenadas, polígonos) de forma integrada com o R.

Vamos agora ler os dados de casos de dengue em São Paulo, obtidos através do TABNET do DATASUS (link):

dengue_sp <- read.csv("https://raw.githubusercontent.com/joaohmorais/dados-mod-temp/
refs/heads/main/csv/dengue_SP.csv"
) head(dengue_sp)
#>   cod_ibge cod_mun6               nome_mun  ano   pop casos
#> 1  3500105   350010             Adamantina 2014 34185    39
#> 2  3500204   350020                 Adolfo 2014  3829     3
#> 3  3500303   350030                  Aguaí 2014 32206   529
#> 4  3500402   350040         Águas da Prata 2014  7565     5
#> 5  3500501   350050       Águas de Lindóia 2014 17623    11
#> 6  3500550   350055 Águas de Santa Bárbara 2014  6171     7
  • O código lê um arquivo .csv do repositório github (um conjunto de dados sobre dengue em São Paulo) e guarda esses dados em um objeto chamado dengue_sp.

  • A função head(dengue_sp), mostra as primeiras linhas da tabela para você visualizar um resumo dos dados carregados.

Temos os dados referente aos municípios, população estimada, o ano de ocorrência, e o número de casos prováveis. A partir dessas informações, podemos calcular a incidência por 10.000 habitantes:

dengue_sp <- dengue_sp %>% 
  mutate(inc = (casos/pop)*10000)

Esse código cria uma coluna chamada inc na tabela dengue_sp, que representa a incidência de dengue como o número de casos (casos) dividido pela população (pop), multiplicado por 10.000.

Em seguida, vamos filtrar somente o ano de 2015 e verificar um resumo da variável incidência:

dengue_2015 <- dengue_sp %>% 
  filter(ano == 2015)
summary(dengue_2015$inc)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    0.00   58.45  135.09  218.71  297.54 1773.28

Esse código faz duas coisas:

  • Ele filtra o conjunto de dados dengue_sp para manter apenas os registros notificados no ano de 2015, criando um objeto chamado dengue_2015.

  • E, por fim, mostra o resumo estatístico, ou seja, mostra estatísticas básicas (mínimo, máximo, média, mediana, etc.) da variável inc (incidência de dengue) do objeto dengue_2015.

Podemos visualizar sua distribuição em um histograma simples:

ggplot(dengue_2015, aes(x=inc)) + 
  geom_histogram(fill = "royalblue", color = "black") + 
  labs(
    x = "Incidência",
    y = "Frequência absoluta",
    title = "Distribuição da Incidência de Dengue por 10.000 \nhabitantes nos municípios de São Paulo, 2015"
  ) +
  theme_minimal() 

Esse código usa a função ggplot2 para criar um histograma que mostra como a incidência de dengue (variável inc) está distribuída nos municípios de São Paulo em 2015.

  • As barras do histograma são definidas como azuis (fill = "royalblue") com contorno preto (color = "black").

  • Os eixos são rotulados como Incidência no eixo dos x e Frequência absoluta n eixo y.

  • O tema visual do gráfico é o minimalista (theme_minimal()), deixando o fundo mais limpo.

Vemos que a distribuição é assimétrica - o que é esperado para uma variável originada por uma contagem - com uma maior concentração de valores mais baixos de incidência, mas com alguns poucos municípios registrando incidências maiores.

Malhas geográficas com geobr

Agora vamos visualizar nossos dados no espaço!

Tradicionalmente, o carregamento de arquivos espaciais para dentro do ambiente R é feito através da leitura de shapefiles (arquivos que incorporam a dimensão espacial aos dados), e usualmente através da função st_read(). Essa função pode ser utilizada da seguinte forma:

malha_exemplo <- st_read("caminho/para/o/shape/malha_br.shp").

Contudo, com a ampla comunidade colaborativa do R, certas operações têm se tornado mais práticas. Hoje, o pacote geobr, desenvolvido pelo IPEA, reune malhas territoriais oficiais de diferentes fontes, como IBGE, CEMADEN e INEP em um repositório único que pode ser acessado diretamente via R, sem necessidade de baixar nenhum shapefile manualmente. Vamos ver como isso funciona:

# se não estiver instalado, rodar:
install.packages("geobr")
library(geobr)

O código está instalando e carregando o pacote geobr no R.

O pacote possui funções específicas para cada nível de agregação espacial, todas começando com read_*. Por exemplo, se quiséssemos uma malha dos estados brasileiros, poderíamos usar read_state(). Para bairros, read_neighborhood(). No caso, iremos recuperar a malha de municípios de São Paulo, e, portanto, utilizaremos read_municipality():

# carregando a malha do estado de São Paulo
malha_sp <- read_municipality("SP", showProgress = F)
head(malha_sp)
#> Simple feature collection with 6 features and 4 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -51.17838 ymin: -23.03648 xmax: -46.54881 ymax: -21.19701
#> Geodetic CRS:  SIRGAS 2000
#>   code_muni              name_muni code_state abbrev_state
#> 1   3500105             Adamantina         35           SP
#> 2   3500204                 Adolfo         35           SP
#> 3   3500303                  Aguaí         35           SP
#> 4   3500402         Águas Da Prata         35           SP
#> 5   3500501       Águas De Lindóia         35           SP
#> 6   3500550 Águas De Santa Bárbara         35           SP
#>                             geom
#> 1 MULTIPOLYGON (((-51.09093 -...
#> 2 MULTIPOLYGON (((-49.69668 -...
#> 3 MULTIPOLYGON (((-47.01254 -...
#> 4 MULTIPOLYGON (((-46.73069 -...
#> 5 MULTIPOLYGON (((-46.635 -22...
#> 6 MULTIPOLYGON (((-49.28903 -...

Esse código carrega o mapa dos municípios do estado de São Paulo (“SP”) e depois exibe as primeiras linhas (informações) desse mapa.

Vemos que os dados estão dispostos de forma similar a uma base de dados tradicional, possuindo agora uma coluna espacial: a coluna geom. Ela contém as informações dos polígonos referentes a cada um dos municípios. É importante notar também o Sistema de Coordenadas (CRS) listado - nesse caso, o SIRGAS 2000.

Podemos verificar rapidamente o contorno dos municípios:

plot(malha_sp$geom)

O código plot(malha_sp$geom) gera um gráfico que desenha a forma geográfica (polígonos ou linhas) que está armazenada no objeto malha_sp, usando a coluna geom, que contém a geometria espacial.

Vemos que em ambos os objetos (dados de dengue e malha espacial) possuímos a coluna referente ao código do município, e, portanto, a usaremos para juntar as duas informações - as informações epidemiológicas e espaciais.

Podemos realizar essa operação através do comando left_join(). Perceba que o código do município é chamado de code_muni na malha espacial e cod_ibge no conjunto de dados de dengue:

malha_sp_dengue <- malha_sp %>% 
  left_join(
    dengue_2015 %>% select(cod_ibge, casos, pop, inc),
    by = c("code_muni" = "cod_ibge")
  )

head(malha_sp_dengue)
#> Simple feature collection with 6 features and 7 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -51.17838 ymin: -23.03648 xmax: -46.54881 ymax: -21.19701
#> Geodetic CRS:  SIRGAS 2000
#>   code_muni              name_muni code_state abbrev_state casos   pop
#> 1   3500105             Adamantina         35           SP  1140 34285
#> 2   3500204                 Adolfo         35           SP   121  3903
#> 3   3500303                  Aguaí         35           SP  2487 32192
#> 4   3500402         Águas Da Prata         35           SP   336  7558
#> 5   3500501       Águas De Lindóia         35           SP   165 17690
#> 6   3500550 Águas De Santa Bárbara         35           SP    16  6313
#>         inc                           geom
#> 1 332.50693 MULTIPOLYGON (((-51.09093 -...
#> 2 310.01793 MULTIPOLYGON (((-49.69668 -...
#> 3 772.55219 MULTIPOLYGON (((-47.01254 -...
#> 4 444.56205 MULTIPOLYGON (((-46.73069 -...
#> 5  93.27304 MULTIPOLYGON (((-46.635 -22...
#> 6  25.34453 MULTIPOLYGON (((-49.28903 -...

O cógido pega a malha dos municípios de SP (malha_sp) e junta (left_join) com uma tabela de casos de dengue de 2015 (dengue_2015), usando como chave primária de ligação o código do município (code_muni em malha_sp e cod_ibge em dengue_2015).

Ele seleciona apenas as colunas cod_ibge, casos, pop (população) e inc (incidência) da tabela dengue_2015 para se vincularem ao objeto malha_sp. No final, o resultado é guardado no objeto malha_sp_dengue.

head(malha_sp_dengue) exibe as primeiras linhas deste novo objeto.

Agora sim! Temos as informações espaciais (coluna geom) e epidemiológicas juntas em um objeto. A partir disso, podemos construir nosso primeiro mapa de incidência de dengue em São Paulo no ano de 2015. Vamos utilizar a função ggplot() que possui uma função específica para visualização de mapas: a geom_sf(). Precisamos especificar duas estéticas (aes()): a geometry, que recebe o nome da coluna que possui as informações espaciais; e fill, que recebe o nome da variável pela qual o mapa será colorido:

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc), linewidth = .1) +
  labs(title = "Incidência de dengue nos municípios de São Paulo, 2015", fill = "Incidência (por 10.000 hab.)") + 
  theme_void()

O código gera um mapa temático mostrando a incidência de dengue nos municípios de São Paulo em 2015.

  • ggplot() inicia o gráfico usando o conjunto de dados malha_sp_dengue, mapeando a geometria dos municípios (geom).

  • geom_sf(aes(fill = inc), linewidth = .1) desenha os municípios coloridos conforme a variável inc (incidência) e determina a largura das linhas do mapa (linewidth = .1).

  • labs() adiciona um título ao mapa e um rótulo para a legenda.

  • theme_void() remove elementos visuais como eixos e grades para destacar apenas o mapa.

Podemos alterar a paleta de cores através dos comandos scale_fill_*, como scale_fill_distiller() e scale_fill_viridis_c() para variáveis contínuas e scale_fill_brewer() e scale_fill_viridis_d() para variáveis discretas ou categorizadas.

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc), linewidth = .1) +
  labs(title = "Incidência de dengue nos municípios de São Paulo, 2015", fill = "Incidência (por 10.000 hab.)") + 
  scale_fill_distiller(palette = "Reds", direction=1) +
  theme_void()

Nesse código utilizado para criar o mapa temático da incidência de dengue nos municípios de São Paulo em 2015 foi utilizado o argumento scale_fill_distiller(palette = "Reds", direction = 1). Ele define o uso da paleta de cores degradê (“Reds”) variando de tons claros a tons escuros de vermelho, seguindo o sentido normal (direction=1), para pintar o mapa segundo os valores de incidência.

Pronto! Temos nosso primeiro mapa de incidência de dengue. Conforme discutido anteriormente, o mapa representa uma taxa contínua, que nem sempre é a melhor opção para visualizar os padrões espaciais - ainda mais quando a distribuição da variável não é simétrica.

Vamos definir, portanto, as quebras - pontos de corte para divisão da incidência no mapa, conforme estratégias que vimos anteriormente.

# quebras pré-definidas (arbitrárias)
quebras_arbitrarias <- c(0, 10, 50, 100, 250, 500, 2000)

malha_sp_dengue <- malha_sp_dengue %>% 
  mutate(
    inc_cat = cut(inc, breaks = quebras_arbitrarias, include.lowest=T)
  )

O código está criando categorias para a variável inc (incidência de dengue) usando faixas de valores pré-definidas (quebras_arbitrarias). Em seguida, adiciona uma nova coluna chamada inc_cat na tabela malha_sp_dengue, indicando a qual intervalo cada valor de inc pertence.

Descrevendo o passo a passo do código:

  • Define as quebras dos intervalos: 0, 10, 50, 100, 250, 500, 2000.

  • Usa a função cut() para categorizar os valores de inc conforme essas faixas.

  • Garante que o menor valor (0) seja incluído no primeiro intervalo (include.lowest = TRUE).

  • Salva essas categorias em uma nova variável chamada inc_cat.

E, ao realizar o gráfico, usamos inc_cat como variável e scale_fill_brewer() para especificar as cores:

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc_cat), linewidth = .1) +
  labs(title = "Incidência de dengue nos municípios de São Paulo, 2015", 
fill = "Incidência (por 10.000 hab.)") + scale_fill_brewer(palette = "Reds", direction=1) + theme_void()

Vamos explorar as outras formas de definir as quebras (fique à vontade para trocá-las no mapa e comparar as diferenças).

# quebras regulares
quebras_regulares <- seq(0,
                         max(dengue_2015$inc), 
                         length.out=7) # length.out define o número de categorias + 1

# quebras quantilicas
quebras_quantilicas <- quantile(
                         dengue_2015$inc, 
                         probs = seq(0, 1, length.out=7))

O código cria dois conjuntos de “quebras” (intervalos) para a variável inc do conjunto de dados dengue_2015:

  1. Quebras regulares (quebras_regulares): Divide o intervalo de inc em 6 categorias de tamanho igual, do valor 0 até o valor máximo de inc.

  2. Quebras quantílicas (quebras_quantilicas): Divide os valores de inc em 6 categorias que contêm aproximadamente o mesmo número de observações, usando os quantis (percentis).

Para as quebras de Jenks (ou quebras naturais), é necessário ter instalada uma implementação do algoritmo de Jenks. Uma opção é através da biblioteca BAMMtools:

install.packages("BAMMtools")
library(BAMMtools)
# quebras naturais ou quebras de Jenks

quebras_jenks <- BAMMtools::getJenksBreaks(
  dengue_sp$inc, k=7
)

quebras_regulares
#> [1]    0.0000  295.5472  591.0944  886.6416 1182.1888 1477.7360 1773.2832
quebras_quantilicas
#>         0%  16.66667%  33.33333%        50%  66.66667%  83.33333%       100% 
#>    0.00000   33.81285   79.98215  135.09475  232.76623  393.97020 1773.28316
quebras_jenks
#> [1]    0.00000   59.01676  174.00383  329.79113  556.57492  933.26643 1773.28316
malha_sp_dengue <- malha_sp_dengue %>% 
  mutate(
    inc_cat2 = cut(inc, breaks = quebras_jenks, include.lowest=T)
  )

O código calcula as quebras naturais (de Jenks) para a variável inc do conjunto dengue_sp, criando 7 classes. Em seguida exibe os três objetos de quebras. E por final o código cria uma nova variável categórica (inc_cat2) na tabela malha_sp_dengue, classificando inc conforme as faixas definidas pelas quebras de Jenks.

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc_cat2) , linewidth = .1) +
  labs(title = "Incidência de dengue nos municípios de São Paulo, 2015", 
) + theme_void()

Pronto! Através da junção de dados epidemiológicos com malhas espaciais construímos nossas primeiras visualizações em mapas com dados de área, alternando paletas e maneiras de definir os pontos de corte.

Vamos voltar agora para o conteúdo teórico do curso!

Autocorrelação espacial

Quando analisamos dados distribuídos espacialmente, há a possibilidade de os dados possuírem autocorrelação espacial. Assim como em séries temporais (veja o curso “Análises de séries temporais em R aplicadas à vigilância em saúde”), a autocorrelação é a presença de correlação, ou dependência, da variável com ela mesmo (no caso temporal, com ela mesma em diferentes pontos no tempo). No contexto espacial, nos referimos à correlação da variável de interesse (como por exemplo, incidência) com ela mesma entre áreas mais próximas.

Ou seja, se há aglomerados de uma certa doença (como vimos nos mapas de dengue), possivelmente os dados possuem uma autocorrelação espacial positiva, pois se um determinado município possui uma alta incidência de dengue, há uma maior probabilidade de seus municípios vizinhos estarem em uma situação parecida.

Há também o caso de autocorrelação inversa, que indica uma correlação inversa entre um município e seus vizinhos. Ou seja, quando seus vizinhos estão com uma incidência baixa, há uma maior probabilidade de o município ter uma incidência alta - esse tipo de autocorrelação pode nos ajudar a identificar áreas destoantes.

Veremos a seguir formas de identificar e mensurar essa autocorrelação presente nos dados, de forma global e local.

Antes de mensurar a presença de dependência espacial nos dados, devemos distinguir esse fenômeno de dependência entre efeitos de primeira e segunda ordem.

Efeitos de primeira ordem - tendência

Quando a dependência espacial sobre os dados se dá devida a efeitos de uma variação de larga escala (ou estrutural), ela é caracterizada como um efeito de primeira ordem. Esses efeitos são observados quando, por exemplo, há tendência de concentração da doença em regiões litorâneas: há uma determinação geográfica subjacente que induz essa tendência, levando a maiores valores de incidência em uma determinada região. Outro caso pode acontecer quando há uma tendência clara de aumento de incidência de uma doença no sentido Norte - Sul.

Vejamos um exemplo de efeito de primeira ordem. Retornando ao nosso conjunto de dados hipotético, vamos supor que uma doença X tenha a seguinte distribuição em um município (suponha população aproximadamente igual entre as áreas):

Figura 31: Exemplo de efeito espacial de primeira ordem: conforme vamos em direção ao leste, os casos da doença X parecem aumentar.

Figura 31: Exemplo de efeito espacial de primeira ordem: conforme vamos em direção ao leste, os casos da doença X parecem aumentar.

Vemos na Figura 31 que a distribuição da doença se dá de forma mais intensa conforme andamos em sentido Leste (ou Sudeste). Pode-se considerar que há, portanto, uma tendência “macro” no processo, em que o mesmo aumenta de intensidade ao avançar em determinado sentido, caracterizando um efeito de primeira ordem. Geralmente, estes efeitos se dão devido a alguma característica ambiental, como maior presença de áreas verdes, maior proximidade com corpos fluviais, maior/menor temperatura, umidade, pressão, entre outras.

Efeitos de segunda ordem - dependência local

Por outro lado, efeitos de segunda ordem (ou efeitos meso, ou efeitos locais) caracterizam uma dependência de uma região sobre seus vizinhos mais próximos. Um exemplo clássico é quando observamos “clusters” de uma doença infecciosa em algumas regiões: se há uma alta de casos em um determinado município, é muito provável que haja também uma alta de casos nos municípios vizinhos, seja pela própria dinâmica infecciosa da doença ou por características ambientais locais propícias para o desenvolvimento de vetores.

Vamos agora nos voltar ao exemplo dos casos de dengue em São Paulo, e verificar a incidência da distribuição de casos em 2020:

dengue_2020 <- dengue_sp %>% 
  filter(ano == 2020)

malha_sp_dengue <- malha_sp %>% 
  left_join(
    dengue_2020 %>% select(cod_ibge, casos, pop, inc),
    by = c("code_muni" = "cod_ibge")
  )

O código filtra os dados de dengue para o ano de 2020 e junta esses dados com a malha dos municípios de São Paulo, criando um objeto chamado malha_sp_dengue que contém informações sobre casos, população e incidência de dengue.

Vamos optar pela classificação das faixas pelas quebras de Jenks:

quebras_jenks <- BAMMtools::getJenksBreaks(
  dengue_2020$inc, k=7
)

malha_sp_dengue <- malha_sp_dengue %>% 
  mutate(
    inc_cat = cut(inc, breaks = quebras_jenks, include.lowest=T)
  )

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc_cat), linewdith = .1) +
  labs(title = "Incidência de dengue nos municípios de São Paulo, em 2020", 
       subtitle = "Quebras de Jenks", 
       fill = "Incidência (por 10.000 hab.)") + 
  scale_fill_brewer(palette = "Oranges", direction=1) +
  theme_void()

Figura 32: Incidência de dengue nos municípios de São Paulo, 2020.

Figura 32: Incidência de dengue nos municípios de São Paulo, 2020.

O código gera um mapa temático da incidência de dengue nos municípios de São Paulo em 2020, usando o método de quebras naturais de Jenks para categorizar os valores.

Passo a passo:

  • Definição das quebras: Calcula 7 intervalos de Jenks (k = 7) para a variável inc (incidência) do conjunto dengue_2020.

  • Criação da variável categórica: No objeto malha_sp_dengue, cria-se uma nova coluna inc_cat, categorizando inc de acordo com as quebras de Jenks.

  • Construção do mapa usando o ggplot2.

Vemos na Figura 32 que há certos “grupos de concentração” da doença no estado. Pode-se listar ao menos quatro regiões em que há um foco de incidência maior entre municípios próximos. Não vemos, no geral, uma tendência clara ao longo de todo o estado como no exemplo anterior. Nesse caso, temos aglomerados da doença nas regiões Oeste e Noroeste do estado, assim como no Sul e Leste também. Esse processo de dependência espacial se aproxima, portanto, de um processo de segunda ordem.

Medidas de proximidade em dados de área

Vimos que as definições de autocorrelação espacial relacionam a variável de interesse entre as áreas próximas ou vizinhas. Portanto, para identificar, testar e quantificar a presença de autocorrelação nos dados de interesse, temos que definir quais áreas são próximas das outras a partir de medidas de proximidade. Essas podem se basear exclusivamente nas fronteiras territoriais, ou em outros fatores como distância (exemplo: distância entre os centroides, ou distância entre suas sedes) ou fluxo (exemplo: fluxo pendular) entre os municípios. Essa escolha se dará a partir dos mecanismos envolvidos no estudo do agravo em questão; há bastante sentido, por exemplo, em considerar o fluxo de deslocamento entre municípios para relacionar espacialmente a distribuição de uma doença como a covid-19. Dessa forma, municípios que tem maior fluxo de pessoas entre si serão mais “próximos” do que municípios geograficamente próximos, mas que não possuem tanto fluxo. Contudo, a disponibilidade de dados sobre esses fatores pode ser um fator limitador, o que nos leva a optar por abordagens mais simples. Vejamos a Figura 33. Neste caso, calculamos a proximidade com base na distância entre os centroides das áreas. Quanto mais perto duas áreas são, mais forte (pois mais próxima) é sua ligação (representada pelas linhas azuis).

Figura 33: Exemplo de cálculo de proximidade entre as áreas do exemplo hipotético.

Figura 33: Exemplo de cálculo de proximidade entre as áreas do exemplo hipotético.

Essas distâncias serão incorporadas a um elemento essencial ao tratar de autocorrelação espacial: a matriz de vizinhança, que veremos a seguir.

Matriz de vizinhança

Uma matriz de vizinhança é uma forma que temos de relacionar as áreas entre si, indicando quais possuem conexões (vizinhança) com as outras e quais os pesos dessas conexões. Geralmente temos uma matriz n×n, onde temos todos os municípios nas linhas e os mesmos municípios nas colunas. Cada célula então, de uma linha l com uma coluna c indica se o município l e c tem relação de vizinhança - e qual o peso dessa relação. Os pesos podem ser fixos (como 0 para uma relação de não vizinhança e 1 para uma relação de vizinhança) ou dependerem de algum outro fator, como a distância entre os municípios.

Vamos considerar o seguinte exemplo, de um recorte de alguns municípios do estado de São Paulo (Figura 34).

Figura 34: Recorte de 7 municípios do estado de São Paulo.

Figura 34: Recorte de 7 municípios do estado de São Paulo.

Uma possível representação das relações entre os municípios exibidos e sua vizinhança seria:

Município Arujá Biritiba-Mirim Guararema Itaquaquecetuba Mogi Das Cruzes Salesópolis Suzano
Arujá 0 0 0 1 1 0 0
Biritiba-Mirim 0 0 1 0 1 1 0
Guararema 0 1 0 0 1 1 0
Itaquaquecetuba 1 0 0 0 1 0 1
Mogi Das Cruzes 1 1 1 1 0 0 1
Salesópolis 0 1 1 0 0 0 0
Suzano 0 0 0 1 1 0 0

Vemos que a primeira linha, que representa o município de Arujá, possui valor 1 nas colunas Itaquaquecetuba e Mogi das Cruzes, pois são os municípios com quem faz fronteira. Mogi das Cruzes, por sua vez, possui valor 1 para quase todos os municípios listados, com exceção de Salesópolis. A matriz exibida é chamada de matriz de vizinhança por contiguidade binária; pois define se um município é vizinho de outro se compartilham fronteira, e caso exista, o peso é atribuído é 1. Caso contrário, o peso é 0. Este tipo de matriz de vizinhança é chamado de vizinhança por contiguidade.

As estratégias mais comuns para definição da matriz de vizinhança são:

  • Por contiguidade: Como no exemplo acima, a célula l,c assumirá valor 1 se o município l e o município c compartilharem fronteira e 0 caso contrário;

  • Por k vizinhos mais próximos: a célula l,c será 1 se c estiver entre os k municípios mais próximos de l. Pode-se considerar, por exemplo, a distância entre as sedes ou centroides dos municípios;

  • Por distância: a célula l,c será 1 se c estiver até a uma distância d pré-definida de l. Caso contrário, será 0.

Vamos implementar cada uma dessas estratégias em nosso mapa do estado de São Paulo.

Por contiguidade

O pacote spdep traz diversas operações envolvendo correlação e dependência espacial implementadas, e, portanto, o utilizaremos para definir nossas matrizes de vizinhança. Caso não esteja instalado, precisaremos rodar install.packages("spdep").

# se não estiver instalado, rodar:
install.packages("spdep")
library(spdep)

A library(spdep) no R é usada para análise de estatística espacial. Ela fornece funções para criar estruturas de vizinhança, calcular pesos espaciais e aplicar testes e modelos que mensuram a autocorrelação espacial.

A partir da malha de municípios de São Paulo, vamos definir nossa matriz de vizinhança. O pacote spdep possui diferentes funções de transformação nesse sentido, onde podemos partir de objetos espaciais como o que temos (malha_sp) e transformá-lo em objetos como de relação entre vizinhos. Uma dessas funções é a poly2nb (polígono (poly) para (to - 2) vizinhança (neighbours - nb)). Ela automaticamente define uma matriz de contiguidade a partir da presença de fronteiras na malha espacial em questão. Vamos testá-la e criar o objeto viz_sp, que utilizaremos mais adiante.

viz_sp <- poly2nb(malha_sp)
viz_sp
#> Neighbour list object:
#> Number of regions: 645 
#> Number of nonzero links: 3530 
#> Percentage nonzero weights: 0.8485067 
#> Average number of links: 5.472868 
#> 1 region with no links:
#> 233
#> 2 disjoint connected subgraphs

A função poly2nb(malha_sp) cria uma lista de vizinhança com base no polígono malha_sp.

  • O objeto malha_sp deve ser um objeto sf ou SpatialPolygons.

  • A função poly2nb() é do pacote spdep.

  • Ela define vizinhos com base no compartilhamento de bordas (isto é, dois polígonos são vizinhos se eles dividem um lado).

  • O objeto viz_sp representa o objeto resultante, que é da classe nb (neighbour list). Cada elemento da lista contém os índices dos polígonos vizinhos de cada polígono de malha_sp.

  • viz_sp impresso no console vai mostrar para cada polígono, quais são os seus vizinhos.

Ótimo! Temos, então, um objeto do tipo “lista de vizinhos” (Neighbour list). Ao chamar o objeto, temos um resumo dessa matriz gerada: são 645 municípios, 3.526 links (ligações), tendo cada região em média 5.47 links (ou seja, 5 vizinhos).

Somos informados também, que uma região (de id 233) não possui links, ou seja, não possui vizinhos. Isso geralmente é um problema quando falamos de análise de dependência espacial. Para testar nossas hipóteses de autocorrelação espacial e utilizar outras técnicas relacionadas, é ideal que todas as regiões tenham ao menos um vizinho. Pode acontecer que, em caso de linhas ou regiões desconexas, ao gerar a lista de vizinhos automaticamente a partir da malha, alguma dessas regiões fique desconectada das demais. Mas lidaremos com isso mais à frente.

Vamos tentar visualizar nossa vizinhança obtida:

# obtendo as coordenadas dos centroides da malha

centroides_sp <- st_centroid(malha_sp)
centroides_coordenadas <- st_coordinates(centroides_sp)

# a partir delas, transformamos nossas relações de vizinhança em linhas:

viz_linhas <- nb2lines(nb = viz_sp, 
                       coords = centroides_coordenadas) %>% 
  st_as_sf()
st_crs(viz_linhas) <- st_crs(malha_sp)

ggplot(malha_sp, aes(geometry=geom)) + 
  geom_sf(fill = "lightblue") + 
  geom_sf(data = viz_linhas, aes(geometry=geometry)) + 
  ggtitle("Matriz de vizinhança por contiguidade, estado de São Paulo.") +
  theme_void()

Figura 35: Matriz de vizinhança por contiguidade dos munícipios do estado de São Paulo.

Figura 35: Matriz de vizinhança por contiguidade dos munícipios do estado de São Paulo.

Esses códigos criam e desenham as relações de vizinhança entre os municípios do estado de São Paulo, representando essas relações com linhas conectando os centroides (pontos centrais) dos polígonos da malha.

  • centroides_sp <- st_centroid(malha_sp): calcula o ponto central (centroide) de cada área (polígono) da malha malha_sp.

  • centroides_coordenadas <- st_coordinates(centroides_sp): Extrai as coordenadas desses centroides para trabalhar com elas depois.

  • viz_linhas <- nb2lines(nb = viz_sp, coords = centroides_coordenadas) %>% st_as_sf(): Converte a lista de vizinhos viz_sp em linhas, usando as coordenadas dos centroides. O resultado é transformado em um objeto sf (simple features).

  • st_crs(viz_linhas) <- st_crs(malha_sp): Define o sistema de referência de coordenadas (CRS) do objeto viz_linhas para ser o mesmo que o da malha original malha_sp.

  • Logo em seguida com a função ggplot, o mapa é plotado primeiro com a malha dos polígonos (malha_sp) com cor azul clara. E depois com as linhas de vizinhança sobre a malha, ligando áreas vizinhas.

Vemos na Figura 35 que todos os municípios que tocam suas fronteiras estão conectados. Os municípios da “borda” do estado tendem a ter menos vizinhos (alguns têm apenas um) e agora conseguimos visualizar o município sem nenhum vizinho. Trata-se de Ilhabela, no Litoral Norte de São Paulo, que justamente por ser uma ilha, não possui fronteira física com nenhum outro município na malha considerada.

Vamos, então, voltar ao nosso objeto de vizinhança viz_sp.

viz_sp
#> Neighbour list object:
#> Number of regions: 645 
#> Number of nonzero links: 3530 
#> Percentage nonzero weights: 0.8485067 
#> Average number of links: 5.472868 
#> 1 region with no links:
#> 233
#> 2 disjoint connected subgraphs

Somos informados que o município que não possui vizinhos é o 233. Vamos verificar se é realmente Ilhabela:

malha_sp$name_muni[233]
#> [1] "Ilhabela"

Pegando o nome do município que está na posição 233 do objeto malha_sp.

Sim, é Ilhabela. Vamos, portanto, conectá-la manualmente a um outro município. Como a chegada a Ilhabela é feita através de balsa partindo de São Sebastião, faz sentindo conectá-la a esse município.

Vamos verificar qual o índice que corresponde a São Sebastião:

which(malha_sp$name_muni == "São Sebastião")
#> [1] 567

Procura a posição (o número do índice) onde o nome do município é “São Sebastião” dentro da coluna name_muni do objeto malha_sp.

Agora que sabemos o índice dos dois municípios, podemos modificar nosso objeto viz_sp para conectá-los:

Incluindo São Sebastião como vizinho de Ilhabela:

viz_sp[[233]] <- 567L

Está forçando a região número 233 (que antes não tinha vizinhos) a ter como vizinho o município de número 567.

Note que utilizamos a notação 567L ao invés de somente 567. Isso se deve porque a estrutura de vizinhança aceita somente índices inteiros. No ambiente R, quando trabalhamos com números sem informar explicitamente seu tipo, o ambiente assuma que é um double (número real de precisão dupla):

class(567)
#> [1] "numeric"
class(567L)
#> [1] "integer"
is.double(567)
#> [1] TRUE
is.double(567L)
#> [1] FALSE

O comando viz_sp[[233]] <- 567L altera manualmente a lista de vizinhança espacial viz_sp, atribuindo ao município de índice 233 o número 567 como seu vizinho. A letra L após o número indica que ele é tratado como um inteiro (integer) no R, o que é necessário porque listas de vizinhança esperam índices inteiros. Essa modificação é útil para corrigir problemas de regiões isoladas, garantindo que todas as áreas tenham pelo menos um vizinho para análises espaciais posteriores.

Logo, precisamos ter essa rigorosidade ao alterar os vizinhos senão o pacote lançará um erro. Agora, vamos fazer o inverso (que também é necessário!): incluir Ilhabela como vizinho de São Sebastião:

viz_sp[[567]] <- c(viz_sp[[567]], 233L)

O comando viz_sp[[567]] <- c(viz_sp[[567]], 233L) está atualizando a lista de vizinhança do município 567, adicionando o município 233 como seu novo vizinho. Ou seja, ele pega os vizinhos que 567 já tinha (viz_sp[[567]]) e junta (c(...)) o número 233 (especificado como inteiro com L), garantindo que a relação de vizinhança seja bidirecional: se 233 tem 567 como vizinho, então 567 também passa a ter 233 como vizinho.

E agora, se repetirmos o comando do gráfico realizado anteriormente, vemos a inclusão da ligação entre Ilhabela e São Sebastião (Figura 36):

viz_linhas <- nb2lines(nb = viz_sp, 
                       coords = centroides_coordenadas) %>% 
  st_as_sf()
st_crs(viz_linhas) <- st_crs(malha_sp)

ggplot(malha_sp, aes(geometry=geom)) + 
  geom_sf(fill = "lightblue") + 
  geom_sf(data = viz_linhas, aes(geometry=geometry)) + 
  labs(title="Matriz de vizinhança por contiguidade, estado de São Paulo.", 
       subtitle = "Após inclusão de Ilhabela") +
  theme_void()

Figura 36: Matriz de vizinhança por contiguidade modificada dos munícipios do estado de São Paulo.

Figura 36: Matriz de vizinhança por contiguidade modificada dos munícipios do estado de São Paulo.

O código acima transforma a lista de vizinhança viz_sp em linhas que conectam os centroides das regiões (centroides_coordenadas) usando a função nb2lines, depois converte essas linhas em um objeto espacial (sf) com st_as_sf(). Em seguida, ajusta o sistema de coordenadas das linhas para ser igual ao da malha malha_sp com st_crs(). Por fim, cria um mapa com ggplot2, desenhando as áreas (malha_sp) preenchidas de azul claro e sobrepondo as linhas de vizinhança (viz_linhas), adicionando título, subtítulo e removendo elementos gráficos extras (theme_void()).

Verificando o objeto de vizinhança novamente:

viz_sp
#> Neighbour list object:
#> Number of regions: 645 
#> Number of nonzero links: 3532 
#> Percentage nonzero weights: 0.8489874 
#> Average number of links: 5.475969 
#> 2 disjoint connected subgraphs

Agora sim! Todos os municípios estão conectados, e não temos nenhum município sem link.

Por vizinhos mais próximos

A segunda estratégia possível que vimos é a definição dos vizinhos a partir dos k municípios mais próximos a um município específico. Repare essa estratégia no exemplo, considerando os dois municípios mais próximos (k=2) para os municípios mostrados anteriormente na Figura 33:

Município Arujá Biritiba-Mirim Guararema Itaquaquecetuba Mogi Das Cruzes Salesópolis Suzano
Arujá 0 0 0 1 0 0 1
Biritiba-Mirim 0 0 0 0 1 1 0
Guararema 0 1 0 0 1 0 0
Itaquaquecetuba 1 0 0 0 0 0 1
Mogi Das Cruzes 0 1 0 0 0 0 1
Salesópolis 0 1 1 0 0 0 0
Suzano 0 0 0 1 1 0 0

Vemos que agora nesse caso podemos ter assimetrias: não necessariamente se c é um dos k municípios mais próximos de l, l será um dos vizinhos mais próximos de c. Isso pode ser visualizado no exemplo apresentado visto que Mogi das Cruzes é vizinha de Guararema (pois está entre suas duas cidades mais próximas) mas Guararema não é vizinha de Mogi, que tem outros dois municípios que são mais próximos (Figura 37).

viz_sp2 <- centroides_coordenadas %>% 
  knearneigh(k=3) %>% 
  knn2nb()

viz_linhas2 <- nb2lines(nb = viz_sp2, 
                       coords = centroides_coordenadas) %>% 
  st_as_sf()
st_crs(viz_linhas2) <- st_crs(malha_sp)

ggplot(malha_sp, aes(geometry=geom)) + 
  geom_sf(fill = "lightgreen") + 
  geom_sf(data = viz_linhas2, aes(geometry=geometry)) + 
  ggtitle("Matriz de vizinhança por k = 3 vizinhos mais próximos, estado de São Paulo.") +
  theme_void()

Figura 37: Matriz de vizinhança por k vizinhos mais próximos dos munícipios do estado de São Paulo.

Figura 37: Matriz de vizinhança por k vizinhos mais próximos dos munícipios do estado de São Paulo.

O script cria uma matriz de vizinhança baseada nos 3 vizinhos mais próximos (k=3) a partir dos centroides dos municípios, desenha as linhas de conexão entre eles e depois plota isso sobre o mapa do estado de São Paulo.

Por distância

Uma terceira estratégia que podemos adotar é a definição de vizinhos a partir de uma distância pré-definida. Nesse caso, a célula l,c será 1 se c estiver até a uma distância d pré-definida de l. Para definir essa distância, é importante entender como as distâncias entre as áreas em questão estão distribuídas.

Vamos exibir a matriz de distâncias e sua distribuição para o nosso exemplo de 7 municípios do estado de São Paulo:

Arujá Biritiba-Mirim Guararema Itaquaquecetuba Mogi Das Cruzes Salesópolis Suzano
Arujá 0,0 40,4 27,1 8,5 24,5 53,4 24,7
Biritiba-Mirim 40,4 0,0 22,2 36,7 17,9 18,9 29,6
Guararema 27,1 22,2 0,0 28,5 20,5 27,9 32,7
Itaquaquecetuba 8,5 36,7 28,5 0,0 19,3 52,0 16,5
Mogi Das Cruzes 24,5 17,9 20,5 19,3 0,0 35,2 13,4
Salesópolis 53,4 18,9 27,9 52,0 35,2 0,0 48,0
Suzano 24,7 29,6 32,7 16,5 13,4 48,0 0,0

Geralmente, trabalhamos com conjunto de áreas muito maiores (estado de São Paulo: 645 municípios), o que torna inviável olhar uma matriz de distância entre cada área (uma matriz 645x645 = 416.025 distâncias). Portanto, é comum recorrermos a métricas do tipo resumo, como médias, medianas e quantis da distribuição de distâncias - ou até mesmo a técnicas gráficas como um histograma - para a definição de um ponto de corte (veremos a seguir). Por enquanto, parece razoável estabelecer o limite de d=20km para definição de nossa matriz de vizinhança com base nas distâncias:

Município Arujá Biritiba-Mirim Guararema Itaquaquecetuba Mogi Das Cruzes Salesópolis Suzano
Arujá 0 0 0 1 0 0 0
Biritiba-Mirim 0 0 0 0 1 1 0
Guararema 0 0 0 0 0 0 0
Itaquaquecetuba 1 0 0 0 1 0 1
Mogi Das Cruzes 0 1 0 1 0 0 1
Salesópolis 0 1 0 0 0 0 0
Suzano 0 0 0 1 1 0 0

Agora, temos uma matriz simétrica novamente - pois se l está a uma distância de até 20km de c, então o contrário também é verdadeiro. Vemos também que o número de vizinhos muda: alguns municípios como Mogi das Cruzes possuem 3 vizinhos pois têm três municípios em até um raio de 20km de distância; outros possuem 1, como Arujá; e Guararema neste caso ficou sem vizinhos, pois está mais afastada dos demais.

Vamos ver como essa situação se aplica ao contexto dos 645 municípios do estado de São Paulo. Primeiro, vamos calcular as distâncias entre os centroides dos municípios. Para isso, convertemos a malha para o sistema de coordenadas UTM (31982), que permite termos coordenadas em metros, ao invés de graus:

# convertendo para UTM, calculando as coordenadas dos centroides novamente e por fim, as distancias:

centroides_coordenadas_metros <- malha_sp %>%
  st_transform(crs = 31982) %>% 
  st_centroid() %>% 
  st_coordinates()

distancias_sp <- centroides_coordenadas_metros %>% 
  dist() 

# um objeto da classe dist
class(distancias_sp)
#> [1] "dist"
#> [1] "dist"
# convertendo para km
distancias_sp <- distancias_sp/1000

distancias_sp %>% 
  as.vector() %>% 
  hist(main = "Histograma das distâncias",
       xlab = "Distância (km)",
       ylab = "Frequência",
       col = "royalblue")

Figura 38: Histograma das distâncias entre os municípios de São Paulo.

Figura 38: Histograma das distâncias entre os municípios de São Paulo.

O script começa transformando a malha espacial malha_sp para o sistema de coordenadas UTM (EPSG 31982), que utiliza metros como unidade de medida. Em seguida, calcula o centroide de cada polígono (isto é, o ponto central de cada município) e extrai as suas coordenadas X e Y.

distancias_sp <- centroides_coordenadas_metros %>% dist(): A partir dessas coordenadas, é calculada a matriz de distâncias entre todos os centroides, resultando em um objeto do tipo dist, onde as distâncias estão inicialmente em metros.

Depois, todas as distâncias são convertidas de metros para quilômetros, dividindo os valores por 1.000.

Por fim, o código transforma essas distâncias em um vetor e constrói um histograma, visualizando a distribuição da frequência das distâncias (em quilômetros) entre os municípios da malha.

Vemos na Figura 38 que a distribuição das distâncias entre os municípios de São Paulo se encontra na faixa de 200-250km, podendo chegar até 800km. Vamos observar os quantis dessa distribuição:

quantile(distancias_sp,
         probs = c(0.01, 0.025, 0.05, 0.10, 0.25, 0.5))
#>        1%      2.5%        5%       10%       25%       50% 
#>  28.12763  44.34926  64.25050  94.52491 163.10356 260.96230

O comando calcula o valor das distâncias que se encontram nos percentis 1%, 2,5%, 5%, 10%, 25% e 50% da distribuição de distâncias entre os municípios.

É ideal que não escolhamos um valor muito alto - escolher 60km, por exemplo, implicará que 5% de todas as relações entre municípios serão consideradas relações de vizinhança - o que pode ser muita coisa. Vamos tentar com o ponto de corte d=40km:

viz_sp3 <- centroides_coordenadas_metros %>% 
  dnearneigh(d1=0, d2=40000)

viz_sp3
#> Neighbour list object:
#> Number of regions: 645 
#> Number of nonzero links: 8464 
#> Percentage nonzero weights: 2.034493 
#> Average number of links: 13.12248

O script cria uma lista de vizinhança espacial chamada viz_sp3 baseada na distância entre centroides.

  • o comando centroides_coordenadas_metros: contém as coordenadas X e Y dos centroides de cada município, em metros (no sistema UTM).

  • dnearneigh(d1=0, d2=40000): Cria a vizinhança espacial considerando que d1 = 0 é a distância mínima entre vizinhos é zero metros (ou seja, qualquer distância positiva serve). E d2 = 40000 é a distância máxima para serem considerados vizinhos é 40.000 metros (40 km).

viz_linhas3 <- nb2lines(nb = viz_sp3, 
                       coords = centroides_coordenadas) %>% 
  st_as_sf()
st_crs(viz_linhas3) <- st_crs(malha_sp)

ggplot(malha_sp, aes(geometry=geom)) + 
  geom_sf(fill = "plum1") + 
  geom_sf(data = viz_linhas3, aes(geometry=geometry),
          alpha=0.4) + 
  ggtitle("Matriz de vizinhança por distância (40km), estado de São Paulo.") +
  theme_void()

Figura 39: Matriz de vizinhança por distância dos munícipios do estado de São Paulo.

Figura 39: Matriz de vizinhança por distância dos munícipios do estado de São Paulo.

O script primeiro cria linhas de vizinhança conectando os centroides dos municípios que são considerados vizinhos segundo o objeto viz_sp3. Essas linhas são convertidas para um formato espacial (sf) para que possam ser usadas em visualizações no ggplot2. Em seguida, o sistema de referência espacial dessas linhas (viz_linhas3) é ajustado para ser o mesmo da malha de municípios (malha_sp), garantindo que os dados estejam sobre a mesma base de coordenadas.

Depois, o script monta o mapa: ele desenha os polígonos dos municípios preenchidos com uma cor lilás clara, sobrepõe as linhas de vizinhança com transparência (para dar destaque ao mapa de fundo), adiciona um título explicativo e remove todos os elementos visuais desnecessários (como eixos, grades e bordas) para deixar o mapa mais limpo e focado na informação espacial.

Vemos na Figura 39 todos os municípios possuem vizinhos, mas algumas áreas são bem mais densas do que outras - municípios rurais tendem a ser maiores e mais espaçados entre uns e outros, enquanto os grandes centros parecem possuir municípios menores e mais próximos um dos outros. Isso se confirma quando vemos o número de vizinhos por município (Figura 40):

# extrair numero de vizinhos de cada municipio

num_vizinhos3 <- viz_sp3 %>% 
  lapply(length) %>% 
  unlist()

malha_sp %>% 
  mutate(num_vizinhos = num_vizinhos3) %>% 
  ggplot(aes(geometry=geom)) + 
  geom_sf(aes(fill=num_vizinhos)) + 
  scale_fill_viridis_c(name="Número de vizinhos",
                       breaks = seq(0,30,by=5)) +
  ggtitle("Número de vizinhos de cada município, segundo a estratégia\nde vizinhança por distância (40km)") +
  theme_void()

Figura 40: Número de vizinhos segundo matriz de vizinhança por distância de 40km, munícipios do estado de São Paulo.

Figura 40: Número de vizinhos segundo matriz de vizinhança por distância de 40km, munícipios do estado de São Paulo.

O script calcula quantos vizinhos cada município possui com base na lista de vizinhança viz_sp3. Primeiro, ele aplica a função length em cada elemento da lista para contar o número de vizinhos de cada município. Esses valores são então organizados em um vetor simples (num_vizinhos3).

Depois, o script usa a malha espacial dos municípios (malha_sp), adiciona uma nova coluna chamada num_vizinhos contendo o número de vizinhos de cada município, e cria um mapa. No mapa, os municípios são coloridos de acordo com a quantidade de vizinhos, utilizando a escala de cores contínuas viridis, que é perceptível para pessoas com daltonismo e facilita a interpretação. O título explica que o critério de vizinhança foi baseado em uma distância de até 40 km entre municípios, e o tema theme_void() é usado para deixar o mapa limpo, sem eixos nem grades.

Dada as diferenças de resultados obtidos, vemos que a escolha da definição da matriz de vizinhança é importante e pode impactar os resultados obtidos nas análises subsequentes. Por isso, é importante ter conhecimento do agravo em questão e entender quais fatores são determinantes para caracterização de uma relação entre um município e outro (divisão de fronteiras, fluxo migratório e/ou pendular, distância, entre outros). Os tópicos a seguir utilizam as matrizes de vizinhança definidas para aplicação de técnicas de análise espacial aos dados.

Testando a autocorrelação espacial global

Durante o capítulo, comentamos sobre o efeito da correlação espacial - o pressuposto de que o espaço - ou seja, quais municípios estão próximos uns dos outros - influencia na determinação do agravo. Esse pressuposto é importante, visto que determinadas abordagens estatísticas - como modelos de regressão - possuem a independência como pressuposto, e caso assumamos independência entre os dados quando há correlação espacial entre eles podemos obter estimativas enviesadas. Mesmo quando o agravo não é diretamente ou intuitivamente relacionado ao espaço, outros fatores determinantes podem se distribuir espacialmente na população e gerar uma dependência espacial do processo de interesse. Por isso, é sempre importante verificar a distribuição dos dados de maneira exploratória e testar a existência dessa correlação.

Geralmente testamos para a autocorrelação espacial - ou seja, a correlação da variável com ela mesma ao longo do espaço - utilizando índices globais conhecidos, como o Índice I de Moran ou o Índice C de Geary.

Índice I de Moran

Para aplicar estes testes, precisamos voltar para a matriz de vizinhança que definimos anteriormente - elas definirão “quem está próximo de quem” e impactarão diretamente em nossas estimativas.

A ideia de cálculo do índice é comparar a variância entre as áreas vizinhas com a variância que seria esperada em caso de aleatoriedade espacial (supondo que o espaço não importa no número de casos do agravo). Os resultados podem ser interpretados da seguinte forma:

  • Índices positivos (I>0): quando I é significativamente maior do que 0, há evidências de que há uma autocorrelação espacial positiva, ou seja: áreas próximas umas das outras tendem a ter taxas parecidas. Geralmente é o caso quando há aglomerados de doença - se um município possui uma taxa alta, é muito provável que seus vizinhos também tenham; quando há baixa atividade da doença, espera-se que os vizinhos também apresentem situação semelhante.
  • Índices negativos (I<0): há uma autocorrelação espacial negativa ou inversa, o que indica que áreas próximas tendem a estar em situações distintas. A definição soa de forma contraintuitiva, mas pode indicar situações em que há outliers presentes: municípios com altas taxas rodeados de municípios com atividade baixa da doença.
  • Índice próximo ou igual a zero (I0): indica ausência de autocorrelação espacial. Os dados estão distribuídos aleatoriamente no espaço, ou seja: um município estar próximo de outro não gera influência no valor do número de casos observado em cada um dos municípios.

No exemplo de casos de dengue em 2020 em São Paulo, há indícios de que exista uma autocorrelação espacial na Figura 32 - visto que há bolsões com maior incidência e concentração da doença. Vamos, portanto, testar a hipótese de autocorrelação espacial.

Para isso, precisamos obter uma matriz de pesos a partir da matriz de vizinhança definida. A partir da relação entre os vizinhos, escolhe-se uma estratégia para ponderar cada uma das ligações (relações de vizinhança) identificadas. É comum a estratégia de ponderar de forma padronizada, de forma que a soma de todos os pesos dos vizinhos de um município m seja igual a 1. Ou seja, se o município m possui dois vizinhos, cada aresta de ligação terá peso 0,5, de forma que ao somar os pesos temos 0,5+0,5=1. Outras estratégias existem e podem ser consultadas com o comando ?nb2listw.

# style = "W" corresponde a abordagem informada de somar as arestas e obter 1 
pesos_viz <- viz_sp %>% nb2listw(style = "W") 

O script cria uma estrutura de pesos espaciais a partir da lista de vizinhança viz_sp, usando o estilo “W” (row-standardized).

A função nb2listw() transforma o objeto de vizinhança em um objeto de pesos (listw), onde:

  • Para cada região, os pesos dos vizinhos são normalizados de forma que a soma dos pesos de cada linha (ou seja, de cada região) seja igual a 1.

  • O argumento style = "W" especifica essa normalização, chamada de “row-standardized weights”.

E agora com os pesos definidos, podemos rodar o teste:

moran.test(
  dengue_2020$inc,
  pesos_viz
)
#> 
#>  Moran I test under randomisation
#> 
#> data:  dengue_2020$inc  
#> weights: pesos_viz    
#> 
#> Moran I statistic standard deviate = 18.034, p-value < 2.2e-16
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic       Expectation          Variance 
#>      0.4408150828     -0.0015527950      0.0006017108

O script realiza o teste de Moran para verificar a presença de autocorrelação espacial na incidência bruta de dengue em 2015.

A função moran.test() recebe:

  • dengue_2020$inc: os valores da variável de interesse (neste caso, a incidência de dengue).

  • pesos_viz: a matriz de pesos espaciais, construída anteriormente com a vizinhança padronizada (“W”).

O teste de Moran avalia se municípios com incidência semelhante estão espacialmente agrupados. A função moran.test() calcula:

  • O valor do índice de Moran (I), que mede a autocorrelação espacial da variável dengue_2020$inc (incidência bruta de dengue).

  • A expectativa do índice sob a hipótese de aleatoriedade espacial (sem padrão).

  • O p-valor associado, que testa se a autocorrelação observada é estatisticamente significativa.

  • Outros detalhes, como variância e método de cálculo.

Observa-se que o Índice I estimado é de 18,03, e o p-valor é menor do que 0,05, indicando a existência de uma autocorrelação espacial positiva. Ou seja, o espaço é um fator importante para a explicação da distribuição da dengue em São Paulo no ano de 2020, de forma que municípios próximos tendem a ter situações semelhantes.

Índice C de Geary

Podemos calcular também o Índice C de Geary, que possui um intuito semelhante, mas se baseia nas diferenças absolutas de valores entre o município e seus vizinhos. Seu valor de neutralidade é o valor 1, com valores de C<1 indicando autocorrelação positiva e C>1 autocorrelação negativa:

geary.test(
  dengue_2020$inc,
  pesos_viz
)
#> 
#>  Geary C test under randomisation
#> 
#> data:  dengue_2020$inc 
#> weights: pesos_viz   
#> 
#> Geary C statistic standard deviate = 13.006, p-value < 2.2e-16
#> alternative hypothesis: Expectation greater than statistic
#> sample estimates:
#> Geary C statistic       Expectation          Variance 
#>       0.544092361       1.000000000       0.001228851

O script realiza o teste de Geary para avaliar a autocorrelação espacial da incidência de dengue em 2015.

A função geary.test() recebe:

  • dengue_2020$inc: a variável de interesse (incidência bruta de dengue).

  • pesos_viz: a matriz de pesos espaciais construída anteriormente.

Ao observar a direção do índice estimado (Geary C statistic < 1) e um p-valor muito baixo, chegamos à mesma conclusão que obtivemos com o teste de Moran.

Ou seja, ambos os índices apontam para uma existência de autocorrelação espacial global nos dados, no sentido positivo. Podemos, contudo, estar interessados em mais: saber onde essa correlação é mais ou menos forte, e se ela muda de sentido ao longo do estado. Para isso, utilizaremos índices locais de associação espacial.

Testando a autocorrelação espacial local

Além de um indicador que estima a existência ou ausência de autocorrelação espacial do processo no conjunto de dados como um todo, pode ser interessante averiguar a existência de fenômenos que ocorrem localmente. Mesmo com ausência de autocorrelação espacial global, podemos observar algumas manifestações locais de autocorrelação, que foram mascaradas ao olhar para a região como um todo. Esses indicadores locais permitem uma exploração mais detalhada das áreas envolvidas em possíveis aglomerados da doença (autocorrelação local positiva), ou a identificação de possíveis outliers (autocorrelação local negativa).

Indicadores Locais de Associação Espacial

O índice de Moran Local (que é um LISA - Local Indicators of Spacial Assocation) é uma adaptação do índice I de Moran para um contexto local, calculado para cada área i. Possui interpretação semelhante, mas voltada para cada área, medindo seu grau de associação com suas áreas vizinhas.

lisa_dengue2020 <- localmoran_perm(dengue_2020$inc, pesos_viz) 

head(lisa_dengue2020)
#>            Ii         E.Ii     Var.Ii        Z.Ii Pr(z != E(Ii))
#> 1  1.40116921  0.010362547 0.20960579  3.03784158    0.002382792
#> 2 -0.24463453 -0.010723014 0.04967033 -1.04954983    0.293925134
#> 3  0.47277760  0.006618033 1.05671851  0.45347667    0.650205542
#> 4  0.02506444  0.003231677 0.10332043  0.06792281    0.945847082
#> 5  0.36274715 -0.005694462 0.12166351  1.05630298    0.290829815
#> 6  0.37898181 -0.007048969 0.07344926  1.42438773    0.154334264
#>   Pr(z != E(Ii)) Sim Pr(folded) Sim   Skewness  Kurtosis
#> 1              0.024          0.012  1.0589117 0.5245390
#> 2              0.300          0.150 -1.0392633 1.2729922
#> 3              0.604          0.302  0.8987260 0.7709809
#> 4              0.744          0.372 -1.7630536 3.2617769
#> 5              0.112          0.056 -1.1696210 0.9156516
#> 6              0.008          0.004 -0.8935898 0.6438619

O script calcula o índice de autocorrelação espacial local (LISA) para a incidência de dengue em 2020 usando permutação aleatória.

A função localmoran_perm():

  • Recebe dengue_2020$inc (incidência bruta) e pesos_viz (estrutura de pesos espaciais).

  • Calcula o índice LISA para cada município, que mede a autocorrelação local — ou seja, identifica clusters de valores altos (hotspots), valores baixos (coldspots) e áreas de transição.

  • Utiliza permutações para gerar uma distribuição nula e calcular p-valores para cada município.

O objeto lisa_dengue2020 armazena:

  • O valor do índice local de Moran para cada área.

  • Estatísticas associadas como p-valor, valores esperados sob aleatoriedade, e variância.

O comando head(lisa_dengue2020) exibe as primeiras linhas, mostrando os resultados dos primeiros municípios analisados.

Agora vemos, para cada município, uma estimativa de I (Ii) e seu p-valor Pr(z != E(Ii)). Podemos representar esse resultado visualmente na Figura 41:

malha_sp_dengue$local_I <- lisa_dengue2020[,1]
malha_sp_dengue$local_I_p_valor <- lisa_dengue2020[,5]

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = local_I)) +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0,
                       name = "I local") +
  theme_void() +
  theme(legend.position = "bottom",
        legend.key.width = unit(1.5, "cm"))

Figura 41: Índice de Moran Local da incidência de dengue nos municípios de São Paulo, 2020.

Figura 41: Índice de Moran Local da incidência de dengue nos municípios de São Paulo, 2020.

O script associa os resultados do LISA (Índice de Moran Local) à malha espacial dos municípios de São Paulo.

Primeiramente:

  • A primeira coluna de lisa_dengue2020 (valores do índice local) é atribuída à nova variável local_I em malha_sp_dengue.

  • A quinta coluna de lisa_dengue2020 (p-valor associado) é atribuída à nova variável local_I_p_valor.

Em seguida, o script cria um mapa utilizando ggplot2:

  • Cada município é preenchido com uma cor baseada no valor do índice local (local_I).

  • Utiliza scale_fill_gradient2(), que aplica uma escala de cores: azul para valores negativos, vermelho para valores positivos e branco para valores próximos de zero.

O mapa resultante destaca visualmente áreas com alta ou baixa autocorrelação local na incidência de dengue em 2015.

Vemos que boa parte do estado concentrou valores de I local próximos de zero. Alguns municípios, no entanto, se destacaram por possuírem valores de I bem maiores que zero, o que significa uma forte autocorrelação espacial local.

Uma forma comum de visualizar os resultados de um LISA é através do chamado LISA map. Nele, utilizamos o p-valor obtido para classificar os municípios com autocorrelação espacial local significativa (p < 0,5) ou não significativa (p >= 0,05). Para os significativos, vemos o sinal da associação: se é direta (I local > 0) ou inversa (I local < 0). Por fim, olhamos para a magnitude da variável: se há uma associação espacial local positiva e a taxa no município é alta (acima da média), então há um cluster Alto-alto, pois é uma região onde altos valores costumam se concentrar. Se a taxa é baixa, contudo, trata-se de um cluster Baixo-baixo, pois há uma concentração de municípios com taxas abaixo do esperado. Quando a associação é significativa e negativa, então há os cenários Alto-baixo e Baixo-alto, que ocorre quando o município é um outlier entre seus vizinhos - se ele está com a taxa alta, seus vizinhos tendem a estar baixo e vice-versa.

No resultado do comando, esse agrupamento já é realizado e armazenado no atributo quadr.

quadrantes <- attr(lisa_dengue2020, "quadr")$mean

malha_sp_dengue$quadrante <- case_when(
  quadrantes == "High-High" ~ "Alto-Alto",
  quadrantes == "Low-Low" ~ "Baixo-Baixo",
  quadrantes == "High-Low" ~ "Alto-Baixo",
  quadrantes == "Low-High" ~ "Baixo-Alto"
)

malha_sp_dengue <- malha_sp_dengue %>% 
  mutate(LISA_resultado = case_when(
    local_I_p_valor < 0.05 ~ quadrante,
    local_I_p_valor >= 0.05 ~ "Não significativo"
  ))

ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = LISA_resultado)) +
  scale_fill_manual(
    values = c("#e41a1c", "#4daf4a", "#377eb8", "gray90")
  ) +
  theme_void() +
  theme(legend.position = "bottom",
        legend.key.width = unit(1.5, "cm"))

Figura 42: LISA Map da incidência de dengue nos municípios de São Paulo, 2020.

Figura 42: LISA Map da incidência de dengue nos municípios de São Paulo, 2020.

O script organiza os resultados do teste LISA em quadrantes e gera um mapa temático para interpretar os padrões espaciais da incidência de dengue em 2020.

Primeiramente:

  • A informação sobre o tipo de associação espacial (quadrante) é extraída dos atributos do objeto lisa_dengue2020, identificando se cada município pertence a “High-High”, “Low-Low”, “High-Low” ou “Low-High”.

  • Cria-se uma nova variável quadrante em malha_sp_dengue, traduzindo esses nomes para português.

Em seguida:

  • Uma nova variável LISA_resultado é criada:

    • Municípios com p-valor menor que 0,05 são classificados no quadrante correspondente (“Alto-Alto”, “Baixo-Baixo”, etc.).
    • Municípios com p-valor maior ou igual a 0,05 são classificados como “Não significativo”.

O mapa final é construído com ggplot2:

  • Cada município é preenchido com uma cor de acordo com o seu resultado no LISA:
    • “Alto-Alto” (vermelho) indica áreas de alta incidência cercadas por outras áreas de alta incidência (hotspots).
    • “Baixo-Baixo” (azul) indica áreas de baixa incidência cercadas por outras de baixa incidência (coldspots).
    • “Alto-Baixo” (verde) e “Baixo-Alto” (verde) indicam áreas de potencial outlier espacial.
    • “Não significativo” (cinza) indica ausência de autocorrelação estatisticamente significativa.
  • O layout do mapa é limpo com theme_void(), e a legenda é posicionada abaixo do mapa para melhor organização visual.

Este mapa LISA permite identificar padrões locais de agregação e outliers espaciais na distribuição da dengue.

Com o resultado, temos os agrupamentos dos municípios de acordo com o resultado do LISA map (Figura 42). Vê-se que a maioria dos municípios não apresentou significância estatística; contudo, alguns aglomerados foram identificados - ressalta-se um cluster de baixa incidência no Sul do estado (grupo “Baixo-Baixo”, em azul), e alguns pontos de clusters de alta incidência (grupo “Alto-Alto”, em vermelho), em regiões já identificadas como destoantes nas visualizações anteriores. O LISA map aponta, inclusive, alguns municípios outliers - que estão com valores de baixa incidência, mas rodeados de municípios com taxas mais altas (grupo “Baixo-Alto”, em verde).

Suavização espacial

Ao representar dados espacialmente, frequentemente convém recorrer a técnicas que facilitam a visualização de padrões espaciais. Isso porque a visualização “crua” dos dados pode conter artefatos que dificultam a percepção e entendimentos dos padrões presentes. Como mencionamos anteriormente, municípios com populações pequenas podem ter “taxas infladas” - caso em que as taxas de incidência de uma doença nesses municípios são muito maiores que os demais devido ao baixo denominador (população) - que faz com que poucos casos (muitas vezes, 2 ou 5 casos) gerem uma taxa muito alta de casos por habitante.

As técnicas de suavização utilizam das relações de vizinhança e proximidade para evidenciar padrões e tendências espaciais, facilitando a visualização e entendimento dos dados. Veremos a seguir duas técnicas diferentes: a suavização por Kernel de Área e pelo Estimador Bayesiano Empírico.

Kernel de área

Como vimos no Módulo 2, podemos interpolar uma superfície contínua de densidade a partir de pontos que representam casos de uma doença. Nesse caso, não temos a localização pontual dos casos - apenas a contagem de casos por municípios, ou outra unidade de área. A aplicação da técnica consiste, portanto, na representação dessas áreas como pontos (geralmente se considera o centroide da área ou alguma outra localização de referência, como a sede do município) - e na ponderação desses pontos com a contagem agregada de casos para cada região.

Para essa transformação de dados de área para pontos, vamos utilizar novamente os centroides dos municípios (Figura 43). Além disso, vamos utilizar o pacote patchwork para organizar as imagens de uma forma mais simples. Caso não tenha instalado, siga os comandos abaixo:

# se não estiver instalado, rodar:
install.packages("patchwork")
library(patchwork)
malha_sp_pontos <- malha_sp %>% 
  st_centroid()

g_malha_sp <- malha_sp %>% 
  ggplot(aes(geometry=geom)) + 
  geom_sf(fill = "darkolivegreen1") + 
  ggtitle("Malha dos municípios de São Paulo") +
  theme_void()

g_pontos_sp <- malha_sp_pontos %>% 
  ggplot(aes(geometry=geom)) + 
  geom_sf(color = "darkgreen") + 
  ggtitle("Municípios de São Paulo, transformados\nem pontos (centroides)") +
  theme_void()

g_malha_sp | g_pontos_sp

Figura 43: Divisão territorial e centroides dos municípios do estado de São Paulo.

Figura 43: Divisão territorial e centroides dos municípios do estado de São Paulo.

O script começa calculando o centroide de cada município da malha malha_sp, ou seja, o ponto central de cada polígono, utilizando a função st_centroid() e criando o objeto malha_sp_pontos.

Em seguida, cria dois gráficos separados com o ggplot2:

  • O primeiro (g_malha_sp) desenha a malha completa dos municípios colorida de verde claro e sem elementos de fundo, apenas os limites dos municípios.

  • O segundo (g_pontos_sp) desenha somente os centroides (pontos centrais dos municípios), com os pontos plotados em preto, também sobre um fundo limpo.

Por fim, usando o pacote patchwork, os dois gráficos são colocados lado a lado para comparação: de um lado aparece o mapa tradicional com os limites dos municípios e, do outro, o mapa reduzido a pontos centrais.

Agora, temos nossos pontos para aplicação do método Kernel. Utilizaremos o Kernel por atributo - onde cada ponto (cada município) terá um valor que corresponde ao que desejamos visualizar, neste caso, a taxa de incidência de uma doença.

Para isso, precisamos criar um objeto da classe ppp (point pattern).

library(spatstat)
library(nngeo)

## Contorno do estado

contorno_sp <- malha_sp %>% 
  st_transform(crs = 31982) %>% 
  st_union() %>% 
  st_remove_holes()

# Transformando o contorno para janela (owin)
sp_owin <- contorno_sp %>% as.owin()

sp_ppp <- ppp(
  centroides_coordenadas_metros[,1], # Coordenada em X
  centroides_coordenadas_metros[,2], # Coordenada em Y
  sp_owin # Janela
)

plot(sp_ppp, pch = 19, cex = 0.5, main = "Objeto PPP dos municípios")

Figura 44: Visualização do objeto de classe ppp contendo os centroides dos municípios e o contorno do estado de São Paulo.

Figura 44: Visualização do objeto de classe `ppp` contendo os centroides dos municípios e o contorno do estado de São Paulo.

O script começa carregando as bibliotecas necessárias para manipulação espacial.

Em seguida, ele cria o contorno do estado de São Paulo: para isso, transforma a malha de municípios para o sistema de coordenadas UTM (metros) através da função st_transform(), une todos os municípios em um único polígono (st_union()) e remove eventuais buracos internos (st_remove_holes()). O contorno gerado servirá para definir a área de referência para análises espaciais futuras, como a estimação de densidades por kernel.

O script primeiro converte o contorno do estado de São Paulo, que estava em formato sf, para o formato owin, usado pelo pacote spatstat para representar áreas de estudo em análises de padrões pontuais.

Em seguida, cria um objeto do tipo ppp, que é a estrutura de dados do spatstat para conjuntos de pontos no espaço. Para isso, utiliza as coordenadas X e Y dos centroides dos municípios e define o contorno do estado como a janela espacial (o limite de estudo).

Por fim, o script plota esse objeto ppp, exibindo os municípios como pequenos pontos pretos dentro do contorno do estado de São Paulo.

Agora temos um objeto do tipo ppp - point pattern, tal qual utilizamos no Módulo 2. Agora, vamos visualizar de forma suavizada no espaço os casos de dengue em 2015 no estado de São Paulo. Podemos voltar a utilizar a função density(), para estimação do kernel, mas teremos que nos atentar a dois parâmetros:

  • sigma: Largura de banda - indicará a ordem de suavização aplicada aos dados (maior ou menor).
  • weights: o parâmetro utilizado para o kernel de atributo. Aqui é importante informarmos o atributo que estamos interessados em suavizar, como por exemplo, a incidência de dengue.
  • eps: indicará a resolução (qualidade gráfica) do resultado. Quanto menor, maior resolução.

Vamos testar as seguintes larguras de banda: 5000, 10000, 20000 e 50000.

kernel_den_2015_b5000 <- density(
  sp_ppp,
  5000,
  weights = dengue_2015$inc,
  scalekernel = T,
  diggle = T,
  eps=2000
)

kernel_den_2015_b10000 <- density(
  sp_ppp,
  10000,
  weights = dengue_2015$inc,
  scalekernel = T,
  diggle = T,
  eps=2000
)

kernel_den_2015_b20000 <- density(
  sp_ppp,
  20000,
  weights = dengue_2015$inc,
  scalekernel = T,
  diggle = T,
  eps=2000
)

kernel_den_2015_b50000 <- density(
  sp_ppp,
  50000,
  weights = dengue_2015$inc,
  scalekernel = T,
  diggle = T,
  eps=2000
)

O script realiza a estimação de densidades por Kernel sobre o objeto de pontos sp_ppp, utilizando diferentes larguras de banda (5000, 10000, 20000 e 50000 metros).

Em cada estimação:

  • A função density() do pacote spatstat é utilizada para calcular a densidade espacial.

  • A largura de banda (sigma) define o grau de suavização: valores menores capturam detalhes locais, enquanto valores maiores suavizam mais amplamente.

  • O argumento weights = dengue_2015$inc pondera os pontos pelo número de casos de dengue em 2015 (inc), ou seja, cada ponto influencia a densidade conforme a intensidade de dengue.

  • scalekernel = TRUE ajusta o kernel para garantir que o total de massa (peso) se mantenha consistente.

  • diggle = TRUE aplica uma correção para bordas, seguindo o método de Diggle, melhorando a estimativa nas áreas próximas às extremidades da janela de estudo.

  • eps = 2000 define a resolução da grade para o cálculo, com espaçamento de 2000 metros entre os pontos de avaliação.

Como resultado, são gerados quatro mapas de densidade (kernel_den_2015_b5000, kernel_den_2015_b10000, kernel_den_2015_b20000, kernel_den_2015_b50000), cada um com um nível diferente de suavização espacial.

Agora, vamos gerar gráficos para compará-los lado a lado. Podemos transformar o resultado em uma tibble():

kernel_den_2015_b5000 <- kernel_den_2015_b5000 %>% 
  as_tibble()
head(kernel_den_2015_b5000)
#> # A tibble: 6 × 3
#>         x        y    value
#>     <dbl>    <dbl>    <dbl>
#> 1 284117. 7498067. 6.26e-14
#> 2 284117. 7500065. 1.77e-13
#> 3 286114. 7494070. 3.57e-14
#> 4 286114. 7496068. 1.39e-13
#> 5 286114. 7498067. 4.60e-13
#> 6 286114. 7500065. 1.30e-12

O script transforma o objeto de densidade kernel_den_2015_b5000, originalmente em formato de grade espacial (im do pacote spatstat), em um tibble, que é uma estrutura de tabela mais moderna e amigável para manipulação no R.

Ao aplicar as_tibble(), a matriz de valores da densidade é convertida em uma tabela onde cada linha representa uma célula da grade espacial, contendo as informações de coordenadas e valor estimado da densidade.

O comando head(kernel_den_2015_b5000) exibe as primeiras linhas dessa tabela, permitindo visualizar como os dados de densidade foram organizados após a conversão.

Nessa tibble, vemos a coordenada estimada em X, Y e o valor de densidade estimado pelo kernel. Lembramos que o kernel sempre estima a densidade do evento em cada localização.

g_kernel_5000 <- kernel_den_2015_b5000 %>% 
  ggplot() + 
  geom_tile(aes(x=x, y=y, fill=value)) + 
  scale_fill_viridis_c(option = "H") +
  labs(title = "sigma: 5000", fill = "densidade") +
  coord_fixed() + theme_void()

kernel_den_2015_b10000 <- kernel_den_2015_b10000 %>% 
  as_tibble()

g_kernel_10000 <- kernel_den_2015_b10000 %>% 
  ggplot() + 
  geom_tile(aes(x=x, y=y, fill=value)) + 
  scale_fill_viridis_c(option = "H") +
  labs(title = "sigma: 10000", fill = "densidade") +
  coord_fixed() + theme_void()

kernel_den_2015_b20000 <- kernel_den_2015_b20000 %>% 
  as_tibble()

g_kernel_20000 <- kernel_den_2015_b20000 %>% 
  ggplot() + 
  geom_tile(aes(x=x, y=y, fill=value)) + 
  scale_fill_viridis_c(option = "H") +
  labs(title = "sigma: 20000", fill = "densidade") +
  coord_fixed() + theme_void()

kernel_den_2015_b50000 <- kernel_den_2015_b50000 %>% 
  as_tibble()

g_kernel_50000 <- kernel_den_2015_b50000 %>% 
  ggplot() + 
  geom_tile(aes(x=x, y=y, fill=value)) + 
  scale_fill_viridis_c(option = "H") +
  labs(title = "sigma: 50000", fill = "densidade") +
  coord_fixed() + theme_void()

(g_kernel_5000 | g_kernel_10000) /
  (g_kernel_20000 | g_kernel_50000)

Figura 45: Estimativa de densidade de kernel da incidência de dengue por município no estado de São Paulo, 2015.

Figura 45: Estimativa de densidade de kernel da incidência de dengue por município no estado de São Paulo, 2015.

O script gera mapas de densidade kernel para quatro larguras de banda diferentes (5000, 10000, 20000 e 50000 metros), convertendo as densidades previamente calculadas para tibbles e visualizando-as com o ggplot2.

Para cada largura de banda:

  • O objeto de densidade (kernel_den_2015_b...) é transformado em tibble usando as_tibble(), organizando as informações de posição (x, y) e valor da densidade (value).

  • Um gráfico (g_kernel_...) é criado usando ggplot(), com geom_tile() para desenhar cada célula da grade como um pequeno quadrado colorido proporcional ao valor da densidade.

  • A escala de cores é definida com scale_fill_viridis_c(option = "H"), que usa um gradiente de cores perceptível e amigável para daltônicos.

  • Títulos e rótulos são adicionados com labs(), informando o valor de sigma utilizado na suavização.

  • coord_fixed() mantém a proporção entre os eixos, e theme_void() remove elementos como eixos e grades para deixar o mapa mais limpo.

Por fim, utilizando o operador de composição de gráficos do patchwork, os quatro mapas (g_kernel_5000, g_kernel_10000, g_kernel_20000 e g_kernel_50000) são organizados em uma grade 2x2 para facilitar a comparação visual dos diferentes níveis de suavização.

Vemos na Figura 45 que a utilização de um sigma baixo, como sigma=5000, gerou uma visualização muito concentrada ainda no centroide dos municípios, suavizando muito pouco e dificultando a percepção de padrões espaciais. Com valores de largura de banda maiores, como sigma=10000 ou sigma=20000 vemos uma suavização interessante da incidência, permitindo ver regiões de maior concentração da doença naquele ano. Vemos o destaque principalmente na região Noroeste de São Paulo, marcada em tons vermelho na escala de cores. O sigma=50000 já realizou uma suavização demasiadamente grosseira, permitindo ver a distribuição em termos macro apenas.

Perceba também, que a escala de densidade muda em cada gráfico: com larguras de banda menores, os tons vermelhos representam uma densidade muito mais alta (9e10-6 = 0,000009), e na maior largura de banda os tons vermelhos representam uma densidade menor (1.5e10-6 = 0,0000015, seis vezes menor).

Vamos comparar a visualização que tínhamos anteriormente da incidência bruta com a visualização suavizada por kernel na Figura 46:

# gráfico feito anteriormente
g_tradicional <- ggplot(malha_sp_dengue, aes(geometry=geom)) + 
  geom_sf(aes(fill = inc)) +
  labs(title = "Incidência bruta", fill = "Incidência (por 10.000 hab.)") + 
  scale_fill_distiller(palette = "Reds", direction=1) +
  theme_void() +
  theme(legend.position = "bottom",
        legend.key.width = unit(1.5, "cm"))

g_kernel <- malha_sp %>% 
  st_transform(crs = 31982) %>% 
  ggplot() + 
  geom_tile(
    data = kernel_den_2015_b10000,
    aes(x=x, y=y, fill = value)
  ) +
  geom_sf(aes(geometry=geom), fill="transparent") + 
  scale_fill_viridis_c(option = "H") +
  labs(title = "Estimativa de densidade de kernel, Sigma=10000", fill = "Densidade") +
  theme_void() +
  theme(legend.position = "bottom",
        legend.key.width = unit(1.5, "cm"))

g_tradicional | g_kernel

Figura 46: Comparação da visualização dos valores burtos e da estimativa de densidade de kernel da incidência de dengue nos municípios do estado de São Paulo, 2015.

Figura 46: Comparação da visualização dos valores burtos e da estimativa de densidade de kernel da incidência de dengue nos municípios do estado de São Paulo, 2015.

O script cria dois gráficos diferentes para comparar a representação da dengue nos municípios de São Paulo em 2015 em relação a incidência bruta e a estimativa de densidade de kernel.

O primeiro gráfico (g_tradicional) utiliza a malha de municípios (malha_sp_dengue) e representa a incidência bruta de cada município usando preenchimento de cor (fill) proporcional ao valor de incidência. O gradiente de cores escolhido é a paleta “Reds”, indicando maior incidência com tons de vermelho mais intensos. A legenda é posicionada na parte inferior e o layout é limpo utilizando theme_void().

O segundo gráfico (g_kernel) utiliza a densidade suavizada gerada pelo método de kernel (kernel_den_2015_b10000). Primeiro, a malha de municípios é transformada para o sistema de coordenadas UTM (metros) usando a função st_transform(). Em seguida, o mapa de densidade é desenhado como uma grade de cores (geom_tile()), e sobre ele é sobreposta a borda dos municípios (geom_sf()) com preenchimento transparente, apenas delineando os limites. A escala de cores é baseada na paleta “viridis” e o sigma utilizado na suavização foi 10000 metros.

Finalmente, utilizando o operador do pacote patchwork, os dois gráficos (g_tradicional e g_kernel) são exibidos lado a lado para permitir uma comparação visual entre a representação tradicional por município e a representação contínua por densidade espacial.

Vemos que, por meio do Kernel por atributo, conseguimos obter uma superfície contínua de densidade de incidência, sobre a qual podemos capturar possíveis padrões espaciais existentes nos dados. Na Figura 46, ao olhar a distribuição da densidade de incidência de dengue em 2015 no estado de São Paulo, notamos a concentração de valores no tom vermelho mais escuro (e mais alto) da escala no Noroeste do estado, além de outros dois ou três aglomerados que atingem os tons amarelados e chamam atenção. Contudo, ao olharmos um dado suavizado por Kernel, olhamos para densidade - o que dificulta a interpretação dos resultados obtidos em termos de casos ou incidência. Veremos a seguir um método de suavização que mantém a mesma unidade da variável de interesse, ou seja, conseguimos suavizar os dados e interpretá-los em termos de incidência.

Método Bayesiano Empírico

Ainda com o intuito de suavização de taxas no mapeamento de doenças, o Método Bayesiano Empírico introduz uma média das taxas que é ponderada de acordo com o tamanho da população do município. Essa média pode ser em relação ao valor global (taxa no estado como um todo, por exemplo) ou local (taxa média nos municípios vizinhos). O método local costuma ser mais interessante no mapeamento de doenças ao fornecer uma suavização mais realista - com base no contexto específico do município e seus vizinhos ao invés de uma média generalizada. A fórmula do método é a seguinte:

θ^i=wiri+(1wi)μi

Perceba que se trata de uma operação de média ponderada entre dois valores (ri e μi), onde o peso é dado por wi.

  • θ^i é a taxa suavizada no município i.
  • ri é a taxa bruta de incidência no município i.
  • μi é a taxa média nos vizinhos de i (ou taxa média geral, se for utilizado o método Bayesiano Empírico Global).
  • e wi é o peso a ser considerado nessa operação de média.

O peso wi é calculado com base nas variâncias da taxa global e da taxa no município i. Quanto maior a população, menor é a variância e assim wi é próximo de 1, dando mais peso à taxa bruta. Quanto menor a população, temos uma variância maior e taxas mais instáveis, e assim mais peso é dado para a média dos vizinhos.

No R temos uma função que realiza esses cálculos diretamente - EBLocal(), do pacote spdep. Vamos aplicá-la ao nosso mesmo caso de incidência de dengue em São Paulo em 2015:

# Relembrando a estrutura dos dados
head(dengue_2015)
#>   cod_ibge cod_mun6               nome_mun  ano   pop casos       inc
#> 1  3500105   350010             Adamantina 2015 34285  1140 332.50693
#> 2  3500204   350020                 Adolfo 2015  3903   121 310.01793
#> 3  3500303   350030                  Aguaí 2015 32192  2487 772.55219
#> 4  3500402   350040         Águas da Prata 2015  7558   336 444.56205
#> 5  3500501   350050       Águas de Lindóia 2015 17690   165  93.27304
#> 6  3500550   350055 Águas de Santa Bárbara 2015  6313    16  25.34453

A função EBLocal() recebe três parâmetros: o número de casos em cada área, a população de cada área, e a estrutura de vizinhança, que definimos anteriormente. Ou seja, vamos optar por uma das estratégias experimentadas (vizinhança por conectividade, vizinhos mais próximos, ou por distância).

inc_bayes_empirico <- EBlocal(
  ri = dengue_2015$casos, 
  ni = dengue_2015$pop,
  nb = viz_sp # matriz de vizinhança definida por contiguidade
)

head(inc_bayes_empirico)
#>           raw         est
#> 1 0.033250693 0.033242118
#> 2 0.031001793 0.030968000
#> 3 0.077255219 0.077160864
#> 4 0.044456205 0.044903289
#> 5 0.009327304 0.009712328
#> 6 0.002534453 0.002692929

O script calcula a estimativa de incidência suavizada usando o método Bayesiano Empírico Local.

A função EBlocal() recebe:

  • ri: o número de casos de dengue observados em cada município (dengue_2015$casos),
  • ni: a população exposta em cada município (dengue_2015$pop),
  • nb: a matriz de vizinhança (viz_sp), que define quais municípios são vizinhos.

O Bayesiano Empírico Local ajusta as taxas observadas considerando a informação dos municípios vizinhos, suavizando variações extremas que podem ser causadas por populações pequenas ou números baixos de casos.

O resultado (inc_bayes_empirico) é um vetor com a taxa de incidência suavizada para cada município. O comando head(inc_bayes_empirico) exibe as primeiras linhas do vetor para inspecionar os valores calculados.

Vemos que o objeto retorna as taxas bruta (raw) e suavizada (est). Vamos extrair a segunda coluna (suavizada) e multiplicá-la por 10.000, como fizemos com a incidência:

dengue_2015 <- dengue_2015 %>% 
  mutate(inc_EB = 1e4*inc_bayes_empirico[,2])

O script adiciona uma nova variável inc_EB ao conjunto de dados dengue_2015.

A operação faz o seguinte:

  • Utiliza mutate() para criar uma nova coluna.

  • inc_bayes_empirico[,2] seleciona a segunda coluna do objeto inc_bayes_empirico, que corresponde à incidência suavizada pelo método de Bayes Empírico.

  • Multiplica esses valores por 10.000 (1e4) para expressar a incidência no formato padrão de casos por 10.000 habitantes.

Assim, a variável inc_EB passa a representar a incidência de dengue suavizada, ajustada para facilitar a comparação com a incidência bruta originalmente calculada.

malha_sp_dengue$inc_EB <- dengue_2015$inc_EB

malha_sp_dengue %>% 
  pivot_longer(cols = c(inc, inc_EB), names_to = "metodo", values_to = "inc") %>%
  mutate(metodo = factor(metodo, 
                         levels = c("inc", "inc_EB"),
                         labels = c("Incidência bruta", "Incidência suavizada por EB Local"))) %>% 
  ggplot(aes(geometry=geom)) + 
  geom_sf(aes(fill = inc)) +
  labs(fill = "Incidência (por 10.000 hab.)") + 
  scale_fill_distiller(palette = "Reds", direction=1) +
  facet_wrap(~metodo) +
  theme_void() +
  theme(legend.position = "bottom",
        legend.key.width = unit(1.5, "cm"),
        plot.title = element_text(hjust=0.5, size=16))

Figura 47: Incidência bruta e suavizada de dengue nos municípios de São Paulo, 2015.

Figura 47: Incidência bruta e suavizada de dengue nos municípios de São Paulo, 2015.

O script compara a incidência bruta de dengue com a incidência suavizada por Bayes Empírico Local (EB Local) nos municípios de São Paulo em 2015.

Primeiramente, a variável inc_EB (incidência suavizada) é adicionada ao objeto espacial malha_sp_dengue.

Depois, utilizando pivot_longer(), os dados são reorganizados de formato “largo” para “longo”, empilhando as duas colunas de incidência (inc e inc_EB) em uma só, com uma nova coluna metodo identificando a origem dos valores.

Em seguida, mutate() ajusta os nomes dos métodos para labels mais amigáveis: “Incidência” para os valores brutos e “Incidência corrigida por EB Local” para os valores suavizados.

O gráfico é construído usando ggplot2, onde:

  • Cada município é desenhado com geom_sf(), preenchido por uma escala de cor proporcional à incidência.

  • scale_fill_distiller() aplica uma paleta de vermelhos para indicar diferentes níveis de incidência.

  • facet_wrap(~metodo) cria dois mapas lado a lado: um para a incidência bruta e outro para a incidência suavizada.

  • O layout é limpo com theme_void(), e elementos como título e legenda são personalizados para melhor apresentação.

O resultado permite comparar visualmente o efeito da suavização por Bayes Empírico Local na distribuição da incidência de dengue.

Vemos que outros padrões chamam a atenção ao comparar a taxa bruta e a taxa suavizada na Figura 47: vê-se que houve redução de algumas flutuações principalmente no sul do estado, e alguns municípios tiveram uma elevação na taxa após a suavização. As maiores diferenças entre incidência bruta e suavizada se dão principalmente nos municípios menores - onde a baixa população implica num menor peso para a incidência bruta calculada, levando a uma maior influência das taxas dos vizinhos.

Principais Modelos de Regressão Espacial para Dados de Área

Além das análises exploratórias para os dados espaciais de área, é possível avançar na modelagem estatística dos dados por meio de modelos que considerem a dependência espacial dos dados. Entre os mais utilizados estão os modelos de autocorrelação espacial global, como o SAR (Simultaneous Autoregressive Model) e o CAR (Conditional Autoregressive Model) (CRESSIE, 1993). O modelo SAR incorpora a dependência espacial diretamente no valor da variável dependente (por exemplo, a incidência de uma doença), assumindo que o valor observado em uma região é influenciado pelos valores das regiões vizinhas. Já o modelo CAR considera essa dependência na estrutura do erro ou resíduo do modelo, sendo muito utilizado em abordagens Bayesianas para modelagem de risco em áreas geográficas.

Outra abordagem é a dos modelos com efeitos espaciais locais, como a Regressão Geograficamente Ponderada (GWR) (FOTHERINGHAM et al., 2002). Diferente dos modelos SAR e CAR, que assumem relações espaciais globais, o GWR permite que os coeficientes da regressão variem ao longo do espaço, ajustando uma equação específica para cada localidade. Isso possibilita identificar como os efeitos de uma variável explicativa sobre a resposta mudam de acordo com a localização, sendo especialmente útil em contextos em que há forte heterogeneidade espacial.

Algumas referências para aprofundamento no tema:

  • CRESSIE, N. (1993). Statistics for Spatial Data (Revised edition). New York: Wiley.

  • FOTHERINGHAM, A. S., BRUNSDON, C., & CHARLTON, M. (2002). Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. Wiley.

Considerações finais

As técnicas exploratórias e os modelos de área são ferramentas cruciais na análise de dados espaciais em epidemiologia, permitindo trabalhar com dados de saúde agregados por unidades geográficas, como municípios, bairros, estados ou distritos sanitários. Sua utilidade é imensa e pragmática, pois alinham a técnicas de visualização e análise estatística à realidade da gestão e vigilância em saúde pública.

A maioria dos dados epidemiológicos da vigilância tais como casos de doenças, óbitos, cobertura vacinal, etc… é disponibilizada publicamente em formato agregado. Isso ocorre tanto por razões operacionais dos sistemas de saúde quanto para proteger a confidencialidade dos indivíduos (conforme a Lei nº 13.709/2018; Lei Geral de Proteção de Dados – LGPD).

A criação dos mapas temáticos que evidenciam o risco de uma doença na região de estudo ajuda a visualizar e identificar padrões espaciais claros — como aglomerados de alto risco (clusters) ou gradientes de risco entre regiões, padrões esse que não seriam visíveis em uma simples tabela de dados.

O uso de modelos espaciais é fundamental e possui diversas aplicações. Algumas delas são:

  1. Permitir a correção de instabilidade em taxas observadas em pequenas áreas com populações pequenas, as quais podem gerar taxas de doença muito instáveis

  2. Identificação de fatores de risco: modelos de regressão, como a regressão de poisson ou binomial negativa, são aplicados para entender quais características da área (variáveis de contexto) estão associadas ao aumento do número de casos. É possível, por exemplo, analisar se municípios com menor cobertura de saneamento básico ou com menor renda média apresentam, de fato, maiores taxas de certas doenças infecciosas.

  3. Ao gerar mapas de risco suavizados e identificar áreas prioritárias, eles permitem que os gestores de saúde pública direcionem recursos de forma mais eficiente, como campanhas de prevenção, fiscalização ou alocação de equipes de saúde, para os locais que mais necessitam.

Essas são apenas alguns exemplos de como os mapas exploratórios e os modelos espaciais podem contribuir com a epidemiologia, vigilância em saúde e tomada de decisões.

Referências

  • BAILEY, T. C. Spatial statistical methods in health. Cadernos de Saúde Pública, Rio de Janeiro, vol. 17, no. 5, p. 1083–1098, Oct. 2001.

  • PEBESMA, E.; BIVAND, R. Spatial Data Science: With Applications in R. 1st ed. New York: Chapman and Hall/CRC, 2023. Available at: https://www.taylorfrancis.com/books/9780429459016. Accessed on: 19 May 2025.

  • PFEIFFER, D. (Ed.). Spatial analysis in epidemiology. Oxford; New York: Oxford University Press, 2008.

  • WALLER, L. A.; GOTWAY, C. A. Applied spatial statistics for public health data. Hoboken, N.J: John Wiley & Sons, 2004(Wiley series in probability and statistics).